如何判断转录组测序是否为链特异性建库

一般而言，我们让公司测序都是知道链特异性建库的方法，除非你觉得他们给的结果不对，或者说你是从公共数据库中获取的转录组数据，并不确定采取的是什么建库方法。

关于链特异性测序的详细介绍，见链特异性测序那点事, 本文只是介绍如何使用RseQC的 infer_experiment.py 命令进行推断而已。

需要注意的是infer_experiment.py 的运行需要一个记录基因注释信息的bed文件，也就是没有注释文件就执行不了。

首先是安装这个RseqQC, 以及准备bedops，可以直接创建一个新的conda环境.

1
2
3
mamba create -n RseqQC -c bioconda -c conda-forge rseqc bedops -y
# 启动环境进行后续分析
conda activate RseqQC

我们需要先将给定的gff文件转成bed文件，假设你的gff文件名是 input.gff

1
gff2bed < input.gff > input.bed12

接着，准备你需要分析样本的BAM文件，如果你只有fastq，那么需要用STAR或者HISAT2进行比对，假设输出结果为 output.bam。

最后就可以使用 infer_experiment.py进行分析

1
2
3
4
5
6
7
infer_experiment.py --input-file output.bam -r input.bed12

# 输出信息如下
This is PairEnd Data
Fraction of reads failed to determine: 0.0160
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0178
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9663

结果表明，read基本都是”1+-,1-+,2++,2–”，也就是最常用的基于dUTP的链特异性建库。

如果是非特异性建库，那么结果大概都是50%-50%。

1
2
3
4
This is PairEnd Data
Fraction of reads failed to determine: 0.0172
Fraction of reads explained by "1++,1--,2+-,2-+": 0.4903
Fraction of reads explained by "1+-,1-+,2++,2--": 0.4925

后续的参数设置，就见链特异性测序那点事