fasta 是一种 基于文本 用于表示 核酸序列 或 多肽序列 的格式。其中核酸或氨基酸均以单个字母来表示,且允许在序列前添加序列名及注释。
特征:2部分-- id行 和 序列行 。
>id行以“>”开头, 后跟序列名称&序列描述。有时候会包含注释信息
>序列行一个字母表示一个 碱基/氨基酸 (A、T、C、G、N (N表示不知道是什么)/20种常见氨基酸)。序列中允许空格,换行,空行,直到下一个“>”,表示该序列结束。
高通量测序(如Illumina NovaSeq等测序平台)得到的原始图像数据文件,经碱基识别(Base Calling)分析转化为原始测序序列(Sequenced Reads),我们称之为Raw Data或Raw Reads,结果以FASTQ(简称为fq)文件格式存储,其中包含测序序列(Reads)的 序列信息 以及其对应的 测序质量信息 。测序样品中真实数据随机截取结果如下图:
特征: 每4行代表一个reads信息
fastq格式是由fasta (记录id和序列) 和QUAL (记录id和碱基质量) 合并而来。fastq文件第三行往往是个+,其实就是和第一行一样都是id。
第四行碱基质量值
碱基质量值(Quality Score或Q-score)是碱基识别(Base Calling)出错的概率的整数映射。通常使用的碱基质量值Q公式[1]为: Q=-10 * log10P 。其中P为碱基识别出错的概率。下表给出了碱基质量值与碱基识别出错的概率的对应关系。
碱基质量值越高表明碱基识别越可靠,准确度越高。比如,对于碱基质量值为Q20的碱基识别,100个碱基中有1个会识别出错,以此类推。
碱基质量值+33(前32个不是单个值),查表找到对应ASCII码
fastq与fasta文件转换
GFF,全称为Generic Feature Format,主要用来描述 基因的结构与功能信息 ,对基因组进行注释。记录序列中转录起始位点、基因、外显子、内含子等组成元件在染色体中的位置信息。现在用得比较多的是第3版,即gff3。gff是一个三级嵌套结构。格式文件为文本文件,分为9列,以TAB分开。控制符使用RFC 3986 Percent-Encoding 编码。比如:%20 代表着ASCII的空格。
gff文件一共有9列:
第九列的详解
GTF全称为gene transfer format,主要是用来对基因进行注释。现在用得比较多的是第2版,即gtf2。gtf文件也是分为9列,前八个字段与GFF相同(有一些小的差别),重点在第九列的不同。
两种文件差异比较:
bam文件和sam文件内容其实是一样的,只是bam是二进制的压缩文件,占内存空间更小。需要通过特定的软件来进行查看。(sam文件可以直接使用 less -S 查看;bam文件使用 samtools view -h xxx.bam | less -S 查看)
SAM(The Sequence Alignment / Map format)格式,即序列比对文件的格式,详细介绍文档: http://samtools.github.io/hts-specs/SAMv1.pdf
SAM文件由两部分组成,头部区和主体区,都以tab分列。
头部区 :以’@'开始,体现了比对的一些总体信息。比如比对的SAM格式版本,比对的参考序列,比对使用的软件等。
主体区 :比对结果,每一个比对结果是一行,有11个主列和一个可选列。
头部区:
@HD VN:1.0 SO:unsorted (排序类型)
头部区第一行:VN是格式版本;SO表示比对排序的类型,有unknown(default),unsorted,queryname和coordinate几种。samtools软件在进行行排序后不能自动更新bam文件的SO值,而picard却可以。
@SQ SN:contig1 LN:9401 (序列ID及长度)
参考序列名,这些参考序列决定了比对结果sort的顺序,SN是参考序列名;LN是参考序列长度;每个参考序列为一行。
例如:@SQ SN:NC_000067.6 LN:195471971
@RG ID:sample01 (样品基本信息)
Read Group。1个sample的测序结果为1个Read Group;该sample可以有多个library的测序结果,可以利用bwa mem -R 加上去这些信息。
例如:@RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq
ID:样品的ID号 SM:样品名 LB:文库名 PU:测序以 PL:测序平台
这些信息可以在形成sam文件时加入,ID是必须要有的后面是否添加看分析要求
@PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7 (比对所使用的软件及版本)
例如:@PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa sampe -a 400 -f ZX1.sam -r @RG ID:ZX1_ID SM:ZX1 LB:PE400 PU:Illumina PL:Miseq …/0_Reference/Reference_Sequence.fa ZX_HQ_clean_R1.fq.sai ZX_HQ_clean_R2.fq.sai …/2_HQData/ZX_HQ_clean_R1.fq …/2_HQData/ZX_HQ_clean_R2.fq
这里的ID是bwa,PN是bwa,VN是0.7.12-r1039版本。CL可以认为是运行程序@RG是上面RG表示的内容,后面是程序内容,这里的@GR内容是可以自己在运行程序是加入的
主体部分介绍:
主体部分有11个主列和1个可选列
FLAG详解: http://broadinstitute.github.io/picard/explain-flags.html
例如:想要查看FLAG 99是什么意思:samtools flags 99
CIGAR详解
CIGAR string,简要比对信息表达式(Compact Idiosyncratic Gapped AlignmentReport),其以参考序列为基础,使用数字加字母表示比对结果,比如3S6M1P1I4M,前三个碱基被剪切去除了,然后6个比对上了,然后打开了一个缺口,有一个碱基插入,最后是4个比对上了,是按照顺序的,字母的含义如下
sam/bam文件查看
samtools工具: http://www.htslib.org/doc/samtools.html
Samtools常用命令的总结:
https://www.bioinfo-scrounger.com/archives/245/
https://www.cnblogs.com/xiaofeiIDO/p/6805373.html
参考: sam格式文件解读
从fasta/fastq文件中提取子集
seqtk subseq head40.fq a.list
提取fq时需要其文件开头用>: sed -i 's/@/>@/g' head40.fq
seqkit grep -f a.list head40.fq [输出格式没有楼上好]
samtools view WT-1.validpairs.bam | head
提取valid_seqid:
提取原始fq_seqid:
提取补集:
提取 validpairs :
提取 discard:
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)