--- title: 【迁移】二代测序数据处理之数据格式说明 urlname: er-dai-ce-xu-shu-ju-chu-li-zhi-shu-ju-ge-shi-shuo-ming date: 2022-01-17 11:33:15 tags: ['NGS', 'fasta', 'fastq', 'gft', 'bam'] --- ## FASTA(.fa) 储存参考数据集 [从零开始完整学习全基因组测序（WGS）数据分析：第2节 FASTA和FASTQ](https://zhuanlan.zhihu.com/p/28470883) * 基本单元 * 序列所表示的基因名：`>ENSMUSG00000020122ENSMUST00000138518`，后可接空格表示注释前缀 * 具体序列信息：`CCCTCCTATCATGC……GGGCCCACCTGTTCTCTGGT` * 基因名独占一行，序列信息为基因名后一行至下一个 `>` 基因名标记前 ```txt >ENSMUSG00000020122ENSMUST00000138518 CCCTCCTATCATGCTGTCAGTGTATCTCTAAATAGCACTCTCAACCCCCGTGAACTTGGT TATTAAAAACATGCCCAAAGTCTGGGAGCCAGGGCTGCAGGGAAATACCACAGCCTCAGT TCATCAAAACAGTTCATTGCCCAAAATGTTCTCAGCTGCAGCTTTCATGAGGTAACTCCA GGGCCCACCTGTTCTCTGGT ``` * FASTA文件为基本单元的简单罗列 ```txt >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)`，上限为 `~` ```txt @DJB775P1:248:D0MDGACXX:7:1202:12362:49613 TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA + JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA ``` ## GTF(.gtf) 描述基因和转录本的信息 [GTF文件格式简介](https://cloud.tencent.com/developer/article/1625204) * 头部有 `#` 开头的注释行 * 主体为 `\t` 分隔的具有九列的表格，空值用 `.` 填充 * 第一列 `seqid` 代表染色体的ID * 第二列是 `source` 代表基因结构的来源 * 第三列是feature, 代表区间对应的特征类型，如外显子等 * 第四、五列为区间的起止坐标 * 第六列是 `score` * 第七列是 `strand`, 代表正负链的信息, +表示正链，-表示负链，?表示不清楚 * 第八列是 `phase`，当描述的是CDS区间信息时，需要指定翻译时开始的位置，取值范围有0,1,2两种 * 第九列是attributes, 表示属性，键值对间以分号分隔，键值对内以空格分隔 ```txt #!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](http://blog.biochen.com/archives/337) * 使用 hisat2-build 工具从.fa文件建立 ```shell 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文件格式详解](https://www.jianshu.com/p/ff6187c97155) bam文件是sam文件的二进制格式，sam 文件是Sequence Alignment/Map Format的简写，产生于比对之后的数据输出，记录了比对的具体情况。文件中以tab键分割，包括 `Header section` 和 `Alignments 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`或者`*`表示。 ```txt @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](https://broadinstitute.github.io/picard/explain-flags.html) #### 第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%。 ## 相关参数说明 [基因组的那些事儿](https://www.jieandze1314.com/post/cnposts/18/) * 测序深度：30x；每个碱基平均被测次数，相关研究表明5~60x中 30x对于后续分析可以达95%置信度 * 测序策略：PE150；PE双端测序、一条序列正反测两次；150每次测150bp，双端测一条片段共300bp * 350bpcDNA建库：将DNA用超声波随机打断成350bp，加接头，作为测序前的准备工作