fastq、fasta、bed、gtf、gff、sam、bam生

2022-01-14 本文已影响0人笺牒九州的怪咖

生信分析过程中的文件格式：
除了原始测序数据fastq、fasta之外，还有基因组文件fasta格式，基因注释文件gtf格式。
在分析的过程中还会有众多中间文件的生成，如bed、bed12、sam、bam、wig、bigwig、bedgraph等。

文件之间的关系与转换方法：

下面来看一下各种文件的详细介绍：

1.测序数据FASTQ文件

1）文件用途：
测序返回的一般数据格式。通常是压缩文件filename.fq.gz的格式。
2）格式说明：
fastq文件每4行代表一条序列
第一行：记录序列测序时所用仪器以及在测序通道中坐标信息，以@开头；
第二行：测序的序列信息，以ATCGN表示，由于荧光信号干扰无法判断是什么碱基时就用N表示；
第三行：通常一个+;
第四行：与第二行碱基信息一一对应，存储测序碱基的质量值（ASCII字符显示）。
3）查看方式：

zcat查看gzip压缩的文件
head -n 8 显示前8行文件内容（前8行代表2条序列）

zcat filename.fq.gz | head -n 8

@SRR1039521.13952745/1
TTCCTTCCTCCTCTCCCTCCCTCCCTCCTTTCTTTCTTCCTGTGGTTTTTTCCTCTCTTCTTC
+
HIJIIJHGHHIJIIIJJJJJJJJJJJJJJJJJJJJJIIJJFIDHIBGHJIHHHHHHFFFFFFE

计算read数
wc -l: 计算行数
bc -l: 计算器 (-l：浮点运算)
echo "zcat N0378901_1.fq.gz | wc-l / (4*1000000)" | bc -l #为什么除以4，又除以1000000，先计算行数再计算的是million值
zcat trt_N0378901_1.fq.gz | awk'{if(FNR%4==0) base+=length}END{print base/10^9,"G";}' # 测序碱基数计算

2.基因组FASTA文件

2.1 基因组FASTA文件

2.1）文件用途
这类fasta文件用于基因组或者基因的DNA或者蛋白的序列信息存储。

2.2）文件格式
以 > 符号开头，记录了该序列类型和所在基因组位置信息，也称“序列名字行”；

“序列行”：一行或者多行，为序列信息，soft-masked基因组会把所有重复区和低复杂区的序列用小写字母标出基因组，小写字母n标示未知碱基。

“不成文”的小规范：
1）第一部分是序列名字，与>相连。
2）第二部分用空格与序列名字相隔，表示注释信息，可以没有。
列如：>gene_00284728 length=231;type=dna

>1 dna_sm:chromosomechromosome:GRCh38:1:1:248956422:1 REF
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
.....
ttgggctggggcctggccatgtgtatttttttaaatttccactgatgattttgctgcatg
gccggtgttgagaatgactgCGCAAATTTGCCGGATTTCCTTTGCTGTTCCTGCATGTA
TTTAAACGAGATTGCCAGCACCGGGTATCATTCACCATTTTTCTTTTCGTTAACTTGCCG
.....
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

通常要求序列名字行简单为好，而且一般加chr作为开头
给第一行添加chr标签，并去掉其他多余信息
下面的写法复杂了些，是为了避免给已经有chr信息的名字再加一次,帮助无脑操作
sed 's/^> \ ([^chr])/ > chr\1/' Homo_sapiens.GRCh38.dna.primary_assembly.fa |cut -f 1 -d ' ' > GRCh38.fa

2.2 测序FASTA文件

2.1）文件用途
这类fasta文件产生于测序的reads的fasta格式，一般为测序的fastq用软件转换而来。

2.2）文件格式
以 > 符号开头，记录了reads序列号，样本信息等；第二行为测得的reads的碱基信息。

>HWI-ST531R 144:D11RDACXX:4:1101:1212:1946 1:N:0:ATTCCT
ATNATGACTCAAGCGCTTCCTCAGTTTAATGAAGCTAACTTCAATGCTGAGATCGTTGACGACATCGAATGGG

3. 基因组注释文件gff和gtf

gff全称General featureformat，主要是用来注释基因组。
gtf全称Gene transfer format，主要是用来对基因进行注释。
两者均是一个9列的基因信息注释文件，前8列的信息几乎一样，区别在于第9列。

gff文件格式：
GFF文件是以tab键分割的9列组成，以下为每一列的对应信息：
1）seq_id：序列的编号，一般为chr或者scanfold编号；
2）source: 注释的来源，一般为数据库或者注释的机构，如果未知，则用点“.”代替。
3）type: 注释信息的类型，比如Gene、cDNA、mRNA、CDS等;
4）start: 该基因或转录本在参考序列上的起始位置；(从1开始，包含);
5）end: 该基因或转录本在参考序列上的终止位置；(从1开始，包含);
6）score: 得分，数字，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测是的P-values值，.表示为空;
7）strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;
8）phase: 仅对注释类型为“CDS”有效，表示起始编码的位置，有效值为0、12 。(对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸，从0开始，CDS的起始位置，除以3，余数就是这个值，，表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置，取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5’末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外；如果Feature为CDS时，必须指明具体值。)；
9）attributes: 一个包含众多属性的列表，格式为“标签＝值”(tag=value)，以多个键值对组成的注释信息描述，键与值之间用“=”，不同的键值用“；”隔开，一个键可以有多个值，不同值用“,”分割。注意如果描述中包括tab键以及“，= ；”，要用URL转义规则进行转义，如tab键用（空格）代替。键是区分大小写的，以大写字母开头的键是预先定义好的，在后面可能被其他注释信息所调用。

预先定义的键主要包括：
ID：注释信息的编号，在一个GFF文件中必须唯一；
name：注释信息的名称，可以重复；Alias：别名；Parent > >
Indicates：该注释所属的注释，值为注释信息的编号，比如外显子所属的转录组编号，转录组所属的基因的编号。
Parent指明feature所从属的上一级ID，用于将exons聚集成transcript，将transripts聚集成gene，值可以为多个；
Target 指明比对的目标区域，一般用于表明序列的比对结果。格式为 “target_idstart end [strand] ,其中strand是可选的(“+”或”-”),target_id中如果包含空格，则要转换成’ ‘。
Gap：T比对结果的gap信息，和Target一起，用于表明序列的比对结果。Derives_from：Note:备注;Dbxref:数据库索引。

gtf文件格式：
GTF格式大部分与GFF相同，但有两个硬性标准：

feature types是必须注明的;
第9列必须以gene_id以及transcript_id开头。而且GTF文件的第9列同GFF文件不同，虽然同样是标签与值配对的情况，但标签与值之间以空格分开，且每个特征之后都要有分号;(包括最后一个特征);
gene_id “geneA”;transcript_id “geneA.1”;database_id “0012”;modified_by “Damian”;duplicates 0;

gtf文件可通过下面的命令对文件进行加工查看：

gunzip Homo_sapiens.GRCh38.94.gtf.gz -c |grep -v '^#' | sed '/ ^ [^ chr ] / s/^/chr/' |less # grep 匹配查询 -v 输出不匹配的行

gff文件与gtf文件相互转换：
使用Cufflinks里面的工具gffread

gff2gtf
gffread my.gff3 -T -o my.gtf

gtf2gff
gffread merged.gtf -o- > merged.gff3

GTF 文件中提取转录本序列(.fa)：

Cufflink中的gffread
gffread transcripts.gtf –g genome.fa –w transcripts.output.fa # 获取转录本序列
gffread transcripts.gtf –g genome.fa -x cds.output.fa # 获取CDS序列
gffread transcripts.gtf –g genome.fa -y protein.output.fa # 获取蛋白序列

Tophat中的gtf_to_fasta
gtf_to_fasta transcripts.gtf genome.fa out_file

从GTF中提取转录本信息

利用AWK命令:
sed 's/"/\t/g' mm10.gtf | awk -F '[\t|;]' '{ OFS="\t" } $3 == "transcript" {print $1,$4,$5,$10,0,$7} ' |sed 's/ transcript_id "//g'|sed 's/"//g'|less -S >mm10.transcript.txt

4. bed文件

bed文件一般代表区域信息，如表示peak位置的bed文件，表示基因注释的bed12文件。

表示基因注释时，gtf/gff和bed文件的区别

1）gtf/gff文件一行表示一个exon/CDS等子区域，多行联合表示一个gene；bed文件一行表示一个gene；
2）gtf文件中碱基位置定位方式是1-based（即起始的碱基记为1），而bed中碱基定位方式是0-based（即起始的碱基记为0）。

bed文件每一行对应信息

必须包含的3列信息：
1）chrom：染色体名字 (e.g.chr3, chrY, chr2_random或者scaffold10671)。
2）chromStart：基因在染色体或scaffold上的起始位置（0-based）。
3）chromEnd：基因在染色体或scaffold上的终止位置（前闭后开）。

可选的9列信息：
4）name：bed文件的行名。
5）score：本条基因在注释数据集文件中的评分（0-1000），在Genome Browser中会根据不同区段的评分显示对应的阴影强度（评分越高灰度越高）。
6）strand：链的方向+、-或. (.表示不确定链的方向)
7）thickStart：CDS区（编码区）的起始位置，即起始密码子的位置。
8）thickEnd：The endingposition at which the feature is drawn thickly (for example the stop codon ingene displays).
9）itemRgb：RGB颜色值（如：255,0,0），方便在GenomeBrowser中查看。
10）blockCount：bed行中外显子的数目。
11）blockSizes：逗号分割的列，数目与blockCount值对应，每个数表示对应外显子的碱基数。
12）blockStarts：逗号分割的列，数目与blockCount值对应，每个数表示对应外显子的起始位置（数值是相对ChromStart计算的）。

5. sam和bam文件

sam文件全称The Sequencing Alignment/Map Format，是Alignment/Map步骤bwa/STAR/HISAT2等软件对结果的标准输出文件，用于存储reads比对到参考基因组的比对结果，是一个纯文本格式，文件一般较大。为了节省硬盘存储，一般使用其高效压缩的二进制格式bam文件。

利用samtools view的-b参数就能把sam文件转为bam文件。

1）sam文件查看方式
在linux终端直接用less即可进行查看；

2）bam文件查看方式
需要借助samtools view工具进行查看

>samtools view filename.bam | less -S
>samtools view -h filename.bam | less -S

NGS分析中大多数文件都是由header和record两部分组成，加上-h参数后可以将header显示出来，默认是不显示的。

@HD    VN:1.5  SO:coordinate@SQ    SN:chr1 LN:248956422@SQ    SN:chr10        LN:133797422......@SQ    SN:chrKI270392.1        LN:971@SQ    SN:chrKI270394.1        LN:970@RG    ID:BH_H3K27ac_2 LB:BH_H3K27ac_2 SM:BH_H3K27ac_2@PG    ID:bwa  PN:bwa  VN:0.7.15-r1140 CL:bwa mem -M -t 8 -R@RG\tID:BH_H3K27ac_2\tLB:BH_H3K27ac_2\tSM:BH_H3K27ac_2\tPL: /MP@PG    ID:MarkDuplicates      VN:1.138(aa51703435dc6a423013e74e56b0b68405facd79_1439324166)   CL:picard.sam.markduplicates.K00141:244:HVL3NBBXX:8:2119:27235:3145399      chr1    10016  32      115M    =      10016   115     CCCTAACCCTAACCCTAACCCK00141:244:HVL3NBBXX:8:2119:27235:31453147     chr1    10016  32      115M    =      10016   -115   CCCTTACCCTAACCCTAACCC

header内容
@HD：是必须的标准文件头，包含版本信息；
@SQ：参考序列染色体名字和长度信息 (SN:scaffold name; LN: length)；
@RG：重要read group信息，通常包含测序平台，测序文库和样本ID等信息，分析时用于区分不同样本（重测序时用到）；
@PG：生成此文件的操作过程和参数信息（program）。
record内容
每一行就是一条read比对上参考基因组的信息，总共12列，用tab键分割。

1. read名称；
2. 比对信息位flag值；
3. 参考序列染色体编号；
4. 5′端起始位置；
5. MAPQ：mapping quality，描述比对的质量，数字越大，特异性越高；
6. CIGAR字符串，记录插入、删除、错配等信息；
7. 配对read所比对到的染色体，仅双端测序的数据才有；
8. 配对read所比对到的位置，仅双端测序的数据才有；
9. 插入片段的长度，仅双端测序的数据才有；
10. read序列；
11. read质量值；
12. 12列以后的信息都是metadata，程序用标记

sam文件中第二列flag信息很重要:

利用samtools flagstat工具可以查看bam文件中比对的flag信息，并输出比对的统计结果。

samtools flagstat *.bam

flag一共有12个标签，使用16进制数表示，每个标签值是2^(n-1)，其中n<=12，每个值有其对应的唯一解释含义:

Decimal Description of read

1 1 Read paired
2 2 Read mapped in proper pair
3 4 Read unmapped
4 8 Mate unmapped
5 16 Read reverse strand
6 32 Mate reverse strand
7 64 First in pair
8 128 Second in pair
9 256 Not primary alignment
10 512 Read fails platform/vendor quality checks
11 1024 Read is PCR or optical duplicate
12 2048 Supplementary alignment

你会发现随机挑选几个值做加和运算，他们的结果都是唯一的，所以在bam文件中第二列flag的值代表这条序列符合下图所示条件的值的和。

所以根据这个值我们可以判断这条序列是双端测序还是单端测序；如果是双端测序，reads来自左端还是右端。比如65 只能是1和64组成，代表这个序列是双端测序，而且是read1。

每次转换很头疼？别担心，网上有很多解码flag含义的在线工具，如SAM Format（网址：https://www.samformat.info/sam-format-flag）.

6. wig、bigwig和bedgraph文件

上述bam和sam文件可以帮助我们探索reads在参考基因组中的比对情况，导入基因组浏览器查看比对状态和突变信息。而wiggle(简称wig)、bigwig(简写bw)以及bedgraph(简写bdg)只包含区域和区域的覆盖度信息，文件更小，用于可视化更方便，可以导入基因组浏览器（Genome Browser）中进行可视化，以查看reads在参考基因组各个区域的覆盖度并检测测序深度。

wiggle：展示区域的密度或强度信息，如GCpercent, probability scores, and transcriptome data.

variableStep chrom=chr2300701 12.5300702 12.5300703 12.5300704 12.5300705 12.5
fixedStep chrom=chr3 start=400601 step=100span=5112233

bedGraph: bed与wig的结合，更省空间和灵活，展示信息与wig类似。(bedGraph的格式一般有四列，Chr、start、end和value，并且坐标是以0为起始左闭右开)

chromA chromStartA chromEndA dataValueAchromB chromStartB chromEndB dataValueB

bigWig: wig文件的二进制压缩格式，可通过wigToBigWig工具转换

推荐大家阅读UCSC官网对这几个文件的详细解释：

wiggle(WIG)：https://genome.ucsc.edu/goldenPath/help/wiggle.html
bedGraph：https://genome.ucsc.edu/goldenPath/help/bedgraph.html
bigWig：http://genome.ucsc.edu/goldenPath/help/bigWig.html

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------I'm a line ! Thanks! --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

参考：
https://mp.weixin.qq.com/s?__biz=MzI5MTcwNjA4NQ==&mid=2247484166&idx=1&sn=417e155672bd718def86003b16bf0078&scene=21#wechat_redirect
https://mp.weixin.qq.com/s/pTyyKfFGUVG2gWFOc4kP_A