转录组分析_PacBio全长转录组分析2

PacBio全长转录组分析2

我们前面已经推送过一篇Pacbio的分析流程了(PacBio全长转录组分析)，但是后来发现对六倍体小麦来说，最后一步是有瑕疵的。另外，也没有校正的介绍，今天我们就补上。听说最近PacBio平台也升级了，号称准确率可以达到99%。

分析流程参考了官方的介绍(https://github.com/PacificBiosciences/IsoSeq_SA3nUP/wiki)。相关命令的安装参见前面一篇的推送。

早期的数据格式有h5格式，包括3个文件，1.bax.h5，2.bax.h5，3.bax.h5。需要先转换成bam格式，命令如下

bax2bam -o mynewbam mydata.1.bax.h5 mydata.2.bax.h5 mydata.3.bax.h5

生成的文件以subreads.bam结尾。

目前再做的话，一般公司交付的就是这个subreads.bam文件。

#获取HiFi reads
ccs -j 10 mydata.subreads.bam mydata.ccs.bam

# Primer removal and demultiplexing
lima --isoseq --peek-guess mydata.ccs.bam primer.fasta mydata.fl.bam

# got FLNC reads
isoseq3 refine mydata.fl.primer_5p--Nanda_3p.bam primer.fasta mydata.flnc.bam --require-polya

# 六倍体小麦，不需要polish。即：不需要 isoseq3 polish
isoseq3 cluster mydata.flnc.bam mydata.polished.bam --verbose --use-qvs

# 将bam文件转换成fastq或者fasta
bam2fasta -o mydata mydata.flnc.bam
bam2fastq -o mydata mydata.flnc.bam

# 使用二代测序reads校正，软件是fmlrc 。首先将同一材料的二代测序reads合并一个文件，双端reads也没关系。
cat mydata.NGS.fastq | awk 'NR % 4 == 2' | sort -T ./ | tr NT TN | ropebwt2 -LR | tr NT TN | fmlrc-convert -f comp_msbwt.npy

fmlrc -p 4 -k 31 -K 61 comp_msbwt.npy mydata.flnc.fasta mydata.flnc.corrected.fa
# mydata.flnc.corrected.fa中存在冗余

最后生成的fasta文件mydata.flnc.corrected.fa里面就是高质量的转录本序列了，接下来就可以进行下游分析了。可以参考这一篇最新PacBio全长转录组标准分析流程。

最后再次提醒下，六倍小麦中最好不要使用isoseq3 polish。因为会把部分同源基因给弄成一个基因。

如果你的数据是公司分析的，最好提醒下这个地方。