二代测序获得的是bcl格式的原始下机数据，通过bcl2fastq软件 https://support.illumina.com/sequencing/sequencing_software/bcl2fastq-conversion-software.html 可将bcl文件转换成每个样本的fq格式文件，也就是我们常理解的数据拆分。bam文件是由比对软件将质控后的fq格式文件与参考基因组进行比对后的比对信息存储文件。
接下来我们理解下bam文件的内容。参考原文提出的一张经典图片：
上图格式的查看方法为：
samtools的header信息每一行都用‘@’ 符号开头，一般大家不会太关注，但其中的信息对于我们有些生信分析还是很重要的。这里需要重点提一下的是header中的@RG也就是Read group信息，这是在做后续数据分析时专门用于区分不同样本的重要信息。比如测序多条lane获得的bam的合并：如果原来样本的测序深度比较深，一般会按照不同的lane分开比对，最后再合并在一起，那么这个时候你会在这个BAM文件中看到有多个RG，里面记录了不同的lane，甚至测序文库的信息，唯一不变的一定是SM的sample信息，这样合并后才能正确处理。这个合并当然也可以在数据拆分后对rawdata进行cat合并，然后再生成bam文件。
接下来是bam的主体内容record（有时候也叫alignment section，即，比对信息），每 一行代表一条reads，每条reads的信息用tab键进行分隔：
对于每列的解释如下表所示：
森或 比如十进制数据77 = 000001001101 = 1 + 4 + 8 +64，这样就得到了这个FLAG包含的意思：PE read，read比对不上参考序列，它的配对read也同样比不上参考序列，它是read1。
二进制的质量描述见下表：
- MAPQ：比对质量值，这个是大家最为熟悉的比对质量值了。比如说Q30（错配率为0.001），Q20（错配率为0.01），计算公式为：-10logP{错比概率} 。一般结果是这一列的数值是从0到60，且0和60这两个数此颂伍字出现次数最多。
- CIGAR：该标签采用数字和几个字符的组合形象记录樱伍了read比对到参考序列上的细节情况，读起来要比FLAG直观友好许多，只是记录的是不同的信息。比如，一条150bp长的read比对到基因组之后，假如看到它的CIGAR字符串为：33S117M，其意思是说在比对的时候这条read开头的33bp在被跳过了（S），紧接其后的117bp则比对上了参考序列（M）。这里的S代表软跳过（Soft clip），M代表匹配（Match）。N表示可变剪接位置，常见于RNA-seq。H只出现在一条read的前端或末端，但不会出现在中间，S一般会和H成对出现，当有H出现时，一定会有一个与之对应的S出现。CIGAR的标记字符有“MIDNSHP=XB”这10个，分别代表read比对时的不同情况：
sam是短序列比对默认的标准格式，是以TAB为分割符的文本格式。主要应用于测序序列mapping到基因组上的结果表示，另外也可以表示其他的多重比对结果。一般把测序reads比对到参考基因组以后，通常得到的就是sam文件。BAM就是SAM的二进制文件，具有更小的存储空间，并且许多下游分析工具使用的是BAM格式。
第1列 ：fastq的read ID
第2列 ：FLAG（如果某一个数值不是下面的任意值，那么那个数值就是下面这些数里面几个的和）
第3列 ：染色体名称。如果这列是“ * ”，可以认为这条read没有比对上的序列，则这一行的第四，五，八，九 列是“0”，第六，七列与该列是相同的表示方法。
第4列 ：比对的位置，从对应上的染色体第1位开始往后计算。没有比对上的，此处为0。
第5列 ：MAPQ比对质量值。越高说明该read比对到参考基因组上的位置越唯一，0表示在参考基因组有多种定位笑纯姿的可能性。60表示在参考基因组只有这一种定位位置。
第6列 ： M表示匹配、I表示插入、D表示删除、N表示内含子和D类似、S表示替换、H表示剪切。
第7列 : 这条reads第二次比对的位置。=表示参考序列与reads一模一样，*表示没有完全一模一样的参考序列。
第8列 : 该列表示与该reads对应的mate pair reads的比对位置（即mate），若无mate,则为0。
第9列 : 序列模板长度，如果同一个片段都比对上了同一个参考序列，为最左边的碱基位置到最右边的碱基位置（左为正，右为负）。当mate 序列位于本序列上游时该值为负值。不可用时，为0。
第10列 : read的序列。
第11列 : ASCII码格式的序列质量。格式同FASTQ一样。其中碰绝1、10、11合起来就是fq格式文件。
第12列 : 可选的区域。格式类似AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:35T0 YT:Z:UU
线性对齐： 一条read比对到参考序列上，可以存在插入(insert)、缺失(delete)、跳跃(skip)、剪切(clip)，但是方向不变（不能是一部分和正链匹配，另一部分又和负链匹配），sam文件中只占用一行记录。
嵌合比对： 由于一条测序read比对到基因组上时分别比对到两个不同的区域，而这两个区域基本没有接触和重叠。因此它在sam文件中需要 占用多行记录 显示。只有第一个记录称作"representative",其他的都是"supplementary"。RNA-seq中裤液的chimeric read或许可以说明有融合基因存在，但在基因组中一般作为结构变异的证据。
---------------------------------------------------------------------------------------------------------------------------------------------------I`m a line ! Thanks !-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
参考链接： https://www.jianshu.com/p/f53741175b67
全长转录组测序（Isoform-sequencing，Iso-seq）基于PacBio单分子实时测序技术（SMRT cell）,凭借超长读长的优势，建库过程中无需打断RNA分子，直接对反转录的全长cDNA测序，得到从5’末端到3’PolyA尾的高质量全长转录本序列，且目前其CCS模式可以达到超高的准确率，可用来进行转录本鉴定、融合基因、可变剪切、精确地分析转录本的丛昌结构等分析。
通过调取polyA尾的全长转录本序列，经反转录成cDNA之后，经过一定规模的扩增，然后进行cDNA损伤修复、末端修复、接头连接、外切酶处理等过程构建Iso-seq RNA文库，其文库构建过程如图一所示。
构建后的哑铃型文库包含测序接头、引物、barcode以及插入片段，如图二所示：
PacBio根据其文库片段长度分为两种模式测序：
其一为 CLR模式 ，对于较长的插入片段，DNA聚合酶的活性不足以支撑合成完全部插入片段，或者仅能合成完一圈多，得到的polymerase reads去除测序接头即为最长subreads；
第二种是 CCS模式 ，对于较短的插入片段，DNA聚合酶的活性可以支撑合成多圈插入片段，此时去除完接头后即为完整的全长的插入片段，同一ZMW孔可产出多个subreads，对subreads进行相互的校验可以得到一致性序列，即CCS（Circular Consensus Sequencing）序列，其10X的准确率可达99.9%，30X可达99.999%。
完整的插入片段序列（Reads of Insert，ROI）一般具有以下特征：包含5’primer、3’primer，且3’primer前存在polyA序列，即Iso-seq文渗卜扒库结构图所示。
因此，理论上我们需要鉴定有这些特征的CCS即可，但实际上，建库过程中会产生嵌合体等非我们需要的序列，需要去过滤掉，整体的转录本鉴定流程可参考图三所示。具体鉴定过程以及实践如下步骤所示。
SMRT cell测序下机后经 smrtlink server初级处理，会将polymerase reads去除接头低质量序列等，转为subreads序列。
具体的用于后续分析的文件为：
movie.subreads.bam
movie.subreads.bam.pbi
movie.subreadset.xml
通过smrttools的ccs工具将subreads.bam转为ccs.bam，具体命令如下：
~/software/smrttools/smrtcmds/bin/ccs movieX.subreads.bam movieX.ccs.bam --min-rq 0.9 # 还可指定--min-passes以及线程数--num-threads
此过程比较耗费资源与时间，如果资源充足，想快速完成ccs的转换，可以对bam文件进行切割，分开转ccs，最后再合并。以下提供了两种并行转ccs的方式，供参考。
获得CCS序列之后，首先需要去掉文库构建过程中的5’和3’测序引物，如果带有barcode，同时也需要去除barcode序列，具体 *** 作可按如下方式：
~/smrttools/smrtcmds/bin/lima movieX.ccs.bam barcoded_primers.fasta movieX.fl.bam --isoseq --peek-guess
其中primer及barcode的格式如下，标签名称必须以“5p”,“3p“结尾，如果有多个3p barcode序列（即包含多个样本），则会同时按照此barcode序列进行拆分，拆分以及去除完引物之后会得到各自样本的bam文件。
文件名称包含引物序弊盯列标签：
movieX.fl.primer_5p--test1_3p.bam
movieX.fl.primer_5p--test2_3p.bam
转录组文库在构建过程中可能会产生嵌合体，即同一个ZMW中两个转录本嵌合到一起。这种嵌合体的出现主要由以下两种情况产生：
鉴于此，这一步需要做的就是对拆分完且去除完引物的CCS序列，进一步过滤，去除嵌合体序列。
~/software/smrtlink/smrtlink_8.0.0.80529/smrtcmds/bin/isoseq3 refine movieX.fl.primer_5p--test1_3p.bam movieX.flnc.bam --require-polya --num-threads 20
由于一个ZMW孔会产生一个转录本序列，即一个CCS，所以不同的CCS可能会是相同的转录本序列，即存在冗余的情况，因此需要再通过聚类（cluster）的方式，对全长转录本序列进行聚类，得到一致性的转录本序列。
Polish纠错是为了进一步提升转录本中碱基的质量，但是这一过程也是非常耗时，目前smrtlink v8版本及以上可以不必进行Polish，即可获得准确度大于0.99的高质量转录本（high-quality isoforms，HQ），和低质量转录本（low-quality isoforms，LQ）。
以上步骤即可得到高质量的转录本序列，其输出结果有如下一些文件。
后续可用polished.hq.fasta.gz进行比对分析等。
全长转录本的鉴定是Iso-seq分析最重要的一步，鉴定出的转录本的质量也决定了后续分析的质量，高质量的转录本可以对转录本的结构进行精确的分析，当然也取决于后续的比对。
参考资料
https://github.com/PacificBiosciences/IsoSeq
https://www.cnblogs.com/xudongliang/p/7473463.html
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）