这一次是主线任务,不过要简单许多。官方帮助:IsoSeq/isoseq-clustering.md at master · PacificBiosciences/IsoSeq · GitHub全英文,悠着点看。
给的原始数据是老师在外面的公司测的,然后给我做一个简单的加工。
我只处理一个P02TYR19102795-1_r64053_20191003_023517_4_H01,所以只需要一路从subreads.bam走到clustered.bam 就算通关了。
STEP0先建立一个虚拟环境,再激活。安装上所需要的软件(也可以一步一步需要的时候再来安装)
conda create -n rnapython=2
conda activate rna
conda install -y pbccs lima pbcoretools
ccs 获取一致性的序列,要求每一条测序的reads至少有一端含有引物序列。但因为导师给我的时候已经是.bam文件了,所以我们跳过第一步。已经拥有了.ccs文件。
STEP1接下来利用lima,把引物从测序结果中剔除出去。引物序列是一个新建的barcoded_primers.fasta文件,应该是测序公司给你的。
有些同学的引物是把brain和liver的样品混在一起测序的,这样子就会有三个序列,和我的稍有不同(但和帮助里面给的代码一样)。
lima能够基于barcode序列区分不同的样品。如果引物(barcoded_primers.fasta)只有两个的话,不要用isoseq。它正是用来去除三个引物时造成的假阳性信号。
接下来就不为引物这个问题纠结了,帮助上全是按有三个来处理的。我们接着进行下一步:Refine。
STEP2这一步的主要目的是去除poly(A),然后把多联体(Concatemer)给去除掉。需要以.fl.bam 或 .fl.consensusreadset.xml和primers.fasta为输入文件。
isoseq3 refine P02TYR19102795-1_r64053_20191003_023517_4_H01.bam barcoded_primers.fasta .flnc.bam --require-polya
但其实我发现输入.bam文件也可以。不用那个.fl。如果要去除poly(A)的话,加上--require-polya。因为我始终是在处理同一批基因,这里的输出就只有.flnc.bam了,后续看着会简便许多。
STEP3接下来做clustering,聚类。生成的clustering就基本是我们需要的东西了。根据纠错结果的不同,lq是得到的低质量序列,hq是高质量序列。
isoseq3 cluster flnc.bam clustered.bam --verbose --use-qvs
还有一个可选操作,对clustering做一个polish
isoseq3 polish clustered.bam P02TYR19102795-1_r64053_20191003_023517_4_H01.subreadset.xml polished.bam
会得到一大堆polish。这一步是将聚类的转录本,合并成一个完整的一致性序列。
但其实还有可以有附加题:统计每一步的Reads。官方给的帮助里面没有,虽然我也找到了相关的东西,等有时间再来做了(如果它真是主线任务的话)。