2022-01-17-【迁移】二代测序数据处理之数据格式说明.md 10 KB


title: 【迁移】二代测序数据处理之数据格式说明 urlname: er-dai-ce-xu-shu-ju-chu-li-zhi-shu-ju-ge-shi-shuo-ming date: 2022-01-17 11:33:15 index_img: https://api.limour.top/randomImg?d=2022-01-17 11:33:15

tags: ['NGS', 'fasta', 'fastq', 'gft', 'bam']

FASTA(.fa) 储存参考数据集

从零开始完整学习全基因组测序(WGS)数据分析:第2节 FASTA和FASTQ

  • 基本单元
  • 序列所表示的基因名:>ENSMUSG00000020122ENSMUST00000138518,后可接空格表示注释前缀
  • 具体序列信息:CCCTCCTATCATGC……GGGCCCACCTGTTCTCTGGT
  • 基因名独占一行,序列信息为基因名后一行至下一个 > 基因名标记前

    >ENSMUSG00000020122ENSMUST00000138518
    CCCTCCTATCATGCTGTCAGTGTATCTCTAAATAGCACTCTCAACCCCCGTGAACTTGGT
    TATTAAAAACATGCCCAAAGTCTGGGAGCCAGGGCTGCAGGGAAATACCACAGCCTCAGT
    TCATCAAAACAGTTCATTGCCCAAAATGTTCTCAGCTGCAGCTTTCATGAGGTAACTCCA
    GGGCCCACCTGTTCTCTGGT
    
  • FASTA文件为基本单元的简单罗列

    >ENSMUSG00000020122ENSMUST00000138518
    CCCTCCTATCATGCTGTCAGTGTATCTCTAAATAGCACTCTCAACCCCCGTGAACTTGGT
    TATTAAAAACATGCCCAAAGTCTGGGAGCCAGGGCTGCAGGGAAATACCACAGCCTCAGT
    TCATCAAAACAGTTCATTGCCCAAAATGTTCTCAGCTGCAGCTTTCATGAGGTAACTCCA
    GGGCCCACCTGTTCTCTGGT
    >……
    ……
    >ENSMUSG00000020122ENSMUST00000125984
    GAGTCAGGTTGAAGCTGCCCTGAACACTACAGAGAAGAGAGGCCTTGGTGTCCTGTTGTC
    TCCAGAACCCCAATATGTCTTGTGAAGGGCACACAACCCCTCAAAGGGGTGTCACTTCTT
    CTGATCACTTTTGTTACTGTTTACTAACTGATCCTATGAATCACTGTGTCTTCTCAGAGG
    CCGTGAACCACGTCTGCAAT
    >……
    ……
    

FASTQ(.fq) 储存原始测序数据

  • 每四行成为一个独立的单元称之为read;FASTQ文件为read的简单罗列
  • 第一行:以‘@’开头,是这一条read的唯一标识符
  • 第二行:测序read的序列,由A,C,G,T和N这五种字母构成,N代表的是测序时那些无法被识别出来的碱基;
  • 第三行:以‘+’开头,用以兼容旧版格式
  • 第四行:测序read的质量值,Q = -10log(测序错误率),字符=chr(ord('!')+Q),上限为 ~

    @DJB775P1:248:D0MDGACXX:7:1202:12362:49613
    TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
    +
    JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA
    

GTF(.gtf) 描述基因和转录本的信息

GTF文件格式简介

  • 头部有 # 开头的注释行
  • 主体为 \t 分隔的具有九列的表格,空值用 . 填充
  • 第一列 seqid 代表染色体的ID
  • 第二列是 source 代表基因结构的来源
  • 第三列是feature, 代表区间对应的特征类型,如外显子等
  • 第四、五列为区间的起止坐标
  • 第六列是 score
  • 第七列是 strand, 代表正负链的信息, +表示正链,-表示负链,?表示不清楚
  • 第八列是 phase,当描述的是CDS区间信息时,需要指定翻译时开始的位置,取值范围有0,1,2两种
  • 第九列是attributes, 表示属性,键值对间以分号分隔,键值对内以空格分隔

    #!genome-build GRCh38.p12
    #!genome-version GRCh38
    #!genome-date 2013-12
    #!genome-build-accession NCBI:GCA_000001405.27
    #!genebuild-last-updated 2018-01
    1       ensembl_havana  gene    65419   71585   .       +       .       gene_id "ENSG00000186092"; gene_version "6"; gene_name "OR4F5"; gene_source "ensembl_havana"; gene_biotype "protein_coding";
    

IDX(.idx) 基因组比对工具HISAT2的索引文件

RNA-Seq基因组比对工具HISAT2

  • 使用 hisat2-build 工具从.fa文件建立

    export FADIR=/opt/human_grch38/dna
    export FANAME=Homo_sapiens.GRCh38.dna.chromosome
    export IDXDIR=/opt/human_grch38/hisat2_idx
    
    export FILELIST=$FADIR/${FANAME}.1.fa,$FADIR/${FANAME}.2.fa,$FADIR/${FANAME}.3.fa,$FADIR/${FANAME}.4.fa,$FADIR/${FANAME}.5.fa,$FADIR/${FANAME}.6.fa,$FADIR/${FANAME}.7.fa,$FADIR/${FANAME}.8.fa,$FADIR/${FANAME}.9.fa,$FADIR/${FANAME}.10.fa,
    export FILELIST=${FILELIST}$FADIR/${FANAME}.11.fa,$FADIR/${FANAME}.12.fa,$FADIR/${FANAME}.13.fa,$FADIR/${FANAME}.14.fa,$FADIR/${FANAME}.15.fa,$FADIR/${FANAME}.16.fa,$FADIR/${FANAME}.17.fa,$FADIR/${FANAME}.18.fa,$FADIR/${FANAME}.19.fa,$FADIR/${FANAME}.20.fa,
    export FILELIST=${FILELIST}$FADIR/${FANAME}.21.fa,$FADIR/${FANAME}.22.fa,$FADIR/${FANAME}.MT.fa,$FADIR/${FANAME}.X.fa,$FADIR/${FANAME}.Y.fa
    echo **************************************
    echo $FILELIST
    echo **************************************
    hisat2-build -p 8 $FILELIST $IDXDIR/GRCh38.hisat2.idx
    

Sam/Bam(.bam) 记录比对的具体情况

Sam/Bam文件格式详解

bam文件是sam文件的二进制格式,sam 文件是Sequence Alignment/Map Format的简写,产生于比对之后的数据输出,记录了比对的具体情况。文件中以tab键分割,包括 Header sectionAlignments section 两部分:

Header section

该部分全部以“@”开头,提供基本的软件版本,参考序列信息,排序信息等

  • @HD行:这一行中有各种不同的标识
    • 标识“VN”用以说明格式版本
    • 标识“SO”用以说明比对排序的情况,有unknown (default)、unsorted、queryname和coordinate,对于coordinate,排序的主键是Alignments section的第三列“RNAME”,其顺序由@SQ行的“SN”标识的顺序定义,次要排序键是Alignments section的第四列“POS”字段。对于RNAME和POS相等的比对,排列顺序则是任意的
  • @SQ行的“SN”标签是参考序列说明,它的值主要是用于Alignments section的第三列“RNAME”和第七列“MRNM”比对的记录
  • @PG行是使用的程序说明;该行“ID”为程序记录标识符,“PN”为程序名字,“CL”为命令行
  • @CO行是任意的说明信息

Alignments section

该部分包含了11列必需字段,无效或者没有的字段一般用0或者*表示。

@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r002 0 ref 9 30 3S6M1P1I4M * 0 0 AAAAGATAAGGATA *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0;
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 6H5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1;
r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1

第1列:Qname

Read的名字

第2列:FLAG

每一个read的比对情况可以用十进制数字(或者十六进制数字)表示,如果比对情况 有多个,将多个比对情况所代表的十进制数字加和就是这一行的FLAG。 另,以下网站可以通过输入FLAG值,直接找出该FLAG是那些FLAG的加和:Decoding SAM flags

第3列:RNAME

比对上的参考序列的名字,该名字出现在Header section的@SQ行的SN标识中,如果该read没有比对上,也就是说该read在参考序列上没有坐标,那么这一列则用“”表示,那么这一行的POS和CIGAR列也会是“”。

第4列:POS

read比对到的参考序列“RNAME”最左侧的位置坐标,也是CIGAR中第一个比对标识“M”对应的最左侧碱基在参考序列的位置,未比对上的read在参考序列中没有坐标,此列标识为“0”。

第5列:MAPQ

比对的质量值,计算方法为比对错误率的-10*log10的值,一般是四舍五入的整数值,如果是255,说明该比对值无效。

第6列:CIGAR

CIGAR标识符表示read中每个碱基的比对情况,主要有以下标识符:

  • M: read上的碱基与参考序列“RNAME”完全匹配,碱基一一对应,包括了正确匹配与错误匹配
  • I: read上的碱基相对于参考序列“RNAME”有插入现象
  • D: read上的碱基相对于参考序列“RNAME”有删除现象
  • N: read上的碱基相对于参考序列“RNAME”存在连续没有比对上的空缺
  • S: read的开头或者结尾部分没有比对到参考序列"RNAME”上, 但这部分未比对上的连续序列仍保留在sam文件的该read序列中
  • H: read的开头或者结尾部分没有比对到参考序列"RNAME”上, 这部分未比对上的连续序列未保留
  • P: padding (silent deletion from padded reference)
  • =:sequence match 正确匹配
  • X:sequence mismatch 错误匹配

第7列:MRNM

该read的mate read比对上的参考序列的名字,该名字出现在Header section的@SQ行的SN标识中,

  • 如果和该read所在行的第三列“RNAME”一样,则用“=”表示,说明这对read比对到了同一条参考序列上;
  • 如果mate read没有比对上,第七列则用“*”表示;
  • 如果这对read没有比对到同一条参考序列,那么这一列则是mate read所在行第三列的“RNAME”。

第8列:MPOS

该read的mate read比对到的参考序列“RNAME”最左侧的位置坐标,也是mate read CIGAR中第一个比对标识“M”对应的最左侧碱基在参考序列的位置,未比对上的read在参考序列中没有坐标,此列标识为“0”。

第9列:ISIZE

表示pair read完全匹配到同一条参考序列时,两个read之间的长度,可简单理解为测序文库的长度。

第10列:SEQ

存储的序列,没有存储,此列则用“*”标识。该序列的长度一定等于CIGAR标识中“M”,“I”,“S”,“=”,“X”标识的碱基长度之和。

第11列:QUAL

序列的每个碱基对应一个碱基质量字符,每个碱基质量字符对应的ASCII码值减去33(Sanger Phred-33 质量值体系),即为该碱基的测序质量得分(Phred Quality Score)。不同Phred Quality Score代表不同的碱基测序错误率,如Phred Quality Score值为20和30分别表示碱基测序错误率为1%和0.1%。

相关参数说明

基因组的那些事儿

  • 测序深度:30x;每个碱基平均被测次数,相关研究表明5~60x中 30x对于后续分析可以达95%置信度
  • 测序策略:PE150;PE双端测序、一条序列正反测两次;150每次测150bp,双端测一条片段共300bp
  • 350bpcDNA建库:将DNA用超声波随机打断成350bp,加接头,作为测序前的准备工作