第一章初步学习生信分析过程中常见文件（fastq/bed/gtf/sam/bam/wig）的格式以及查看方式

最新推荐文章于 2024-05-07 19:55:49 发布

Laurel*

最新推荐文章于 2024-05-07 19:55:49 发布

阅读量6.2k

点赞数 9

分类专栏：生信爱好者文章标签：生信初体验

生信爱好者专栏收录该内容

0 篇文章 6 订阅

订阅专栏

   本人是从今年的四月份起开始接触生信，刚开始学也是一头雾水，根本不知道如何下手，自己在网上查看别人发的博客和各种资料，看Linux和R语言书籍，再一步一步跟着别人做，走了不少弯路，也学到了不少东西，现在才算是稍微懂了一点点，所以想把自己学到的一点点生信知识分享给更多初步接触生信的同学，尤其是作RNA-seq和CHIP-seq的同学们，同时也希望大家能够解决我的一些困惑，共同学习，共同进步，哈哈。另外，由于刚开始写博客，一方面是为了能够与同行的朋友交流学习，另一方面更多是为了记录自己的一点小经验。刚开始接触生信，我们必须要对其常见的一些文件格式搞清楚，这样才有便于我们更好的理解每一步分析的意义在哪里，每一步分析是否正确，是否符合要求。所以，本章主要是整合了他人的博客，然后自己在ncbi上下载了一组数据进行了一下试用，原创性低，大家适当参考即可。但是，整体跑完后觉得这些代码很好理解也很适合初学者，故特作此分享，有不足之处敬请谅解！

测序数据FASTQ文件

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

第一行以@开头，后面是reads的ID以及其他信息，例如上例中 HWUSI-EAS100R代表Illmina设备名称，6代表flowcell中的第六个lane，73代表第六个lane中的第73个tile，941:1973代表该read在该tile中的x：y坐标信息；#0，若为多样本的混合作为输入样本，则该标志代表样本的编号，用来区分个样本中的reads；/1代表paired end中的前一个read。第二行为read的序列。第三行以“+”开头，跟随者该read的名称（一般于@后面的内容相同），但有时可以省略，但“+”一定不能省。第四行代表reads的质量
Q值得计算
Illumina测序仪是按照荧光信号来判断所测序的碱基是哪一种的，例如红黄蓝绿分别对应ATCG，那么一旦出现一个紫色的信号该怎么判断呢，因此对每个结果都有一个概率的问题。起初sanger中心用Phred quality score来衡量该read中每个碱基的质量，既-10lgP ，其中P代表该碱基被测序错误的概率，如果该碱基测序出错的概率为0.001，则Q应该为30，那么30+33=63，那么63对应的ASCii码为“？”，则在第四行中该碱基对应的质量代表值即为“？”
ASCII码对照表
http://tools.jb51.net/table/ascii
其他常用命令
https://blog.csdn.net/qazplm12_3/article/details/85222665

基因组FASTA文件
类似于fastq文件格式，这里不再重述，详细内容参考https://baike.so.com/doc/6690561-6904467.html
基因组注释文件gff和gtf
gff全称General feature format，主要是用来注释基因组。gtf全称Gene transfer format，主要是用来对基因进行注释。两者均是一个9列的基因信息注释文件，前8列的信息几乎一样，区别在于第9列。具体参考https://mp.weixin.qq.com/s/rZ26i19hiS5ZOqIoqkL1Wg
从ensemble下载的gtf文件前5行一般是以#开头的注释信息，后续分析中用不上需要去除，同时需要给第一列添加chr标签（与基因组序列一致），可通过下面的命令对文件进行加工：

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

表示基因注释时，gtf/gff和bed文件的区别gtf/gff文件一行表示一个exon/CDS等子区域，多行联合表示一个gene；bed文件一行表示一个gene；gtf文件中碱基位置定位方式是1-based，而bed中碱基定位方式是0-based
zcat ./GCF_000001635.26_GRCm38.p6_genomic.gff.gz
4. bed文件
分析过程中的bed文件一般代表区域信息，如表示Peak位置的bed文件，表示基因注释的bed12文件。
必须包含的3列信息：
1）chrom：染色体名字 (e.g. chr3, chrY, chr2_random或者scaffold10671)。
2）chromStart：基因在染色体或scaffold上的起始位置（0-based）。
3）chromEnd：基因在染色体或scaffold上的终止位置（前闭后开）。
可选的9列信息：
name：bed文件的行名。
score：本条基因在注释数据集文件中的评分（0-1000），在Genome Browser中会根据不同区段的评分显示对应的阴影强度（评分越高灰度越高）。
strand：链的方向+、-或. (.表示不确定链的方向)
thickStart：CDS区（编码区）的起始位置，即起始密码子的位置。
thickEnd：The ending position at which the feature is drawn thickly (for example the stop codon in gene displays).
itemRgb：RGB颜色值（如：255,0,0），方便在Genome Browser中查看。
blockCount：bed行中外显子的数目。
blockSizes：逗号分割的列，数目与blockCount值对应，每个数表示对应外显子的碱基数。
blockStarts：逗号分割的列，数目与blockCount值对应，每个数表示对应外显子的起始位置（数值是相对ChromStart计算的）。
less ./SRX091905MC_peaks.bed
5. sam和bam文件
sam文件查看方式
在linux终端直接用less即可进行查看；
bam文件查看方式
需要借助samtools view工具进行查看
samtools view filename.bam | less -S或samtools view -h filename.bam | less -S
NGS分析中大多数文件都是由header和record两部分组成，加上-h参数后可以将header显示出来，默认是不显示的：
@HD：是必须的标准文件头，包含版本信息；@SQ：参考序列染色体名字和长度信息 (SN: scaffold name; LN: length)@RG：重要read group信息，通常包含测序平台，测序文库和样本ID等信息，分析时用于区分不同样本（重测序时用到）；@PG：生成此文件的操作过程和参数信息（program）
record内容

在这里插入图片描述
sam文件中第二列**flag**信息很重要，下面做进一步解释，利用samtools flagstat工具可以查看bam文件中比对的flag信息，并输出比对的统计结果。在这里插入图片描述

samtools flagstat SRX091905BW.bam或 samtools flagstat <in.bam> |<in.sam> | <in.cram>
每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以“PASS + FAILED”格式显示。
flag一共有12个标签，使用16进制数表示，每个标签值是2^(n-1)，其中n>=12（这个没看懂），每个值有其对应的唯一解释含义，具体见下图
在这里插入图片描述

你会发现随机挑选几个值做加和运算，他们的结果都是唯一的，所以在bam文件中第二列flag的值代表这条序列符合下图所示条件的值的和
所以根据这个值我们可以判断这条序列是双端测序还是单端测序；如果是双端测序，reads来自左端还是右端。比如65 只能是1和64组成，代表这个序列是双端测序，而且是read1

samtools flags 的含义https://blog.csdn.net/xcaryyz/article/details/79257327
在这里插入图片描述

Laurel*

关注

9
点赞
踩
53

收藏

觉得还不错? 一键收藏
0
评论
第一章初步学习生信分析过程中常见文件（fastq/bed/gtf/sam/bam/wig）的格式以及查看方式

测序数据FASTQ文件# zcat查看gzip压缩的文件# head -n 8 显示前8行文件内容（前8行代表2条序列）zcat filename.fq.gz | head -n 8 Q值得计算Illumina测序仪是按照荧光信号来判断所测序的碱基是哪一种的，例如红黄蓝绿分别对应ATCG，那么一旦出现一个紫色的信号该怎么判断呢，因此对每个结果都有一个概率的问题。起初sanger中...
复制链接

扫一扫