利用RseQC进行判断，其python小脚本infer_experiment.py可以帮助我们判断链特异性测序。
准备bed文件：
convert2bed
convert2bed是BEDOPS软件包里非常常用的函数，可以把常见的二进制或者文本基因组文件（BAM, GFF, GTF, GVF, PSL, OUT, SAM, VCF，WIG）转换成bed格式。这里的bed格式大部分都是bed6格式。
conda install -c bioconda bedops
convert2bed -i gff < A.gff3 > A.bed
此时为6位bed，还要转换为12位
bed6Tobed12.py
从Ensembl数据库网站下载的gtf或gff文件，利用convertbed可以轻松转换成bed文件，但是不能直接转成bed12格式。笔者利用python3写了一个脚本bed6Tobed12.py，可以实现将bed6文件转换成bed12格式。
脚本根据bed6文件中的第四列ID值，将ID值相同的行抽取出来，根据基因组坐标排序，计算行数及每行相对第一行的起始位置，转换成bed12格式输出。
注意，相同ID的行之间基因组坐标不能相互重合，重合的行可以先利用bedtools merge合并，再转换成bed12格式。
# 下载脚本
wget https://github.com/ustbcaoqi/biotools/archive/master.zip
unzip master.zip
# 进行转换
python3 biotools-master/bed6Tobed12.py file.bed(6)>file_new.bed(12)
得到12位bed文件
安装rseqc（可尝试python2环境安装）：
conda install -c bioconda rseqc
安装完成后发现，其实并没有安装软件，而是多了一些脚本，其中有一个就叫infer_experiment.py。
使用
基本使用格式如下：
Usage: infer_experiment.py [options]
Options:
--version show program's version number and exit
-h, --help show this help message and exit
-i INPUT_FILE, --input-file=INPUT_FILE
Input alignment file in SAM or BAM format
-r REFGENE_BED, --refgene=REFGENE_BED
Reference gene model in bed fomat.
-s SAMPLE_SIZE, --sample-size=SAMPLE_SIZE
Number of reads sampled from SAM/BAM file.
default=200000
-q MAP_QUAL, --mapq=MAP_QUAL
Minimum mapping quality (phred scaled) for an
alignment to be considered as "uniquely mapped".
default=30
部分参数：
（1） -i 比对生成的bam文件（可以不用排序）
（2） -r gtf转bed12文件产生的bed文件。技能——gtf转为bed12
（3） -s 从所有的reads中抽取多少进行统计（默认200k）
（4） -q unique map的mapq阈值
双端实例
infer_experiment.py -i D0-1.bam -r ../bed12.bed -q 20
Reading reference gene model /share/home/ShuiKM/reference/refdata-GRCm39/GRCm39.genes.bed ... Done
Loading SAM/BAM file ... Total 200000 usable reads were sampled
This is PairEnd Data
Fraction of reads failed to determine: 0.0712
Fraction of reads explained by "1++,1--,2+-,2-+": 0.4649
Fraction of reads explained by "1+-,1-+,2++,2--": 0.4640
这个双端测序的数据显然是非链特异性的，详细来看：
1++,1--,2+-,2-+表示read1比对上的链和基因所在的链一样，read2比对上的链和基因所在的链是相反的；
1+-,1-+,2++,2--表示read1比对上的链和基因所在的链是相反的，read2比对上的链和基因所在的链一样。
问题就在于这两种情况的占比是几乎一样的，说明read的链与基因所在的链完全是随机的，显然就不是链特异性测序了。
链特异性
转自作者：Bio_Infor、清昭_QCao
转自链接：https://www.jianshu.com/p/109997345a67
链接：https://www.jianshu.com/p/847801e8bf92
网友评论