从NCBI上下载转录组数据，很多文章方法描述简单，无法判断是否为链特异性数据，导致在mapping和raw reads count时参数不知如何选择。一直弄不懂链特异性和参数设置的同志们可以去看《链特异建库那点事》，讲的非常清楚（虽然小白看完还是一知半解）。所以关于判断转录组数据是否为链特异性，偷懒的小白找到了一个比较省事儿的方法，用RSeQC的infer_experiment.py工具。
官网链接：http://rseqc.sourceforge.net/#infer-experiment-py
官网用法如下：
infer_experiment.py的用法举例
需要输入自己数据的bam文件，这个容易拿到。
但是另一个hg19.refseq.bed12文件，官网给出的参数解释如下：
infer_experiment.py的参数
小白不知道Reference gene model in bed format是啥，呜呜呜～，于是在官网开始找相关信息。目录为Input format一栏显示了bed12文件的举例文件和推荐的使用工具：
input format内容
但是Bedops(Bedopts)的gff2bed工具和RSeQC的举例bed12格式却大不一样。
Bedops(Bedopts)的gff2bed转化成bed文件的结果：
原链接：
https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/gtf2bed.html
gff2bed结果
而RSeQC举例的bed12格式的文件，除了包含常用的chromosome, start, end, name, score, strand等信息外，最后一列包含了多个extron和intron等的位置，用逗号隔开。
原链接：
http://dldcc-web.brc.bcm.edu/lilab/liguow/RSeQC/dat/sample.bed
RSeQC举例的bed12文件
小白通过搜索终于找到了获取bed12文件Reference gene model in bed format的方法，就是使用UCSC的gtfToGenePre工具，小白在上一篇笔记《Reference gene model in bed format》中已经详细讲述，这里只列代码：
#安装gtfToGenePre
conda install -c bioconda ucsc-gtftogenepred
#准备好基因组gtf文件，从gtf转换为GenePred格式
gtfToGenePred -genePredExt -geneNameAsName2 genes.gtf gene.tmp
#从GenePred文件提取信息得到bed文件
awk '{print $2"\t"$4"\t"$5"\t"$1"\t0\t"$3"\t"$6"\t"$7"\t0\t"$8"\t"$9"\t"$10}' gene.tmp > genes_refseq.bed12
拿到bed12文件后，开始试试用infer_experiment.py判断数据是否为链特异性。
infer_experiment.py -r genes_refseq.bed12 -i col.bam
结果：
Reading reference gene model genes_refseq.bed12 ... Done
Loading SAM/BAM file ... Total 200000 usable reads were sampled
This is SingleEnd Data
Fraction of reads failed to determine: 0.0050
Fraction of reads explained by "++,--": 0.4974
Fraction of reads explained by "+-,-+": 0.4976
这表明该数据为单端数据，以illumina standard建库方式为代表的fr-unstranded的非链特异性转录组数据。
再来一个链特异性的双端数据试试：
infer_experiment.py -r genes_refseq.bed12 -i m1_col_1_s.bam
结果：
Reading reference gene model genes_refseq.bed12 ... Done
Loading SAM/BAM file ... Total 200000 usable reads were sampled
This is PairEnd Data
Fraction of reads failed to determine: 0.0057
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0049
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9894
这表明该数据为双端数据，以dUTP建库方式为代表的RF (fr-firststrand)的链特异性转录组数据。
结果的详细解读可以去官网或参考链接：
http://rseqc.sourceforge.net/#infer-experiment-py
https://www.jianshu.com/p/4987dce4d165
欢迎关注和讨论哦，小白们一起学习生信～
网友评论