因为需要统计不同RNA_type对应的clean reads长度,首先需要获得不同RNA对应的gtf文件,转化为bed文件,再从我们真实数据的bam文件中提取出相应的bam,根据提取的bam从包含有clean reads的fastq文件中提取出clean reads,最后统计clean reads的长度,作图。
PS:不知道有什么巧妙的办法能够快速统计的,欢迎大佬多多交流。
有一种方法俺找到咧,samtools stats,里面LR、FLR、LRL都可以根据bam文件统计,very方便。
1. 生成bed文件
1.1 生成gtf文件
- 获取RNA type和gene name的对应关系,这里是根据ensembl网站上下载的数据先得到gene name对应的RNA类型。这里因为工作需要大概分成了十三种,可以得十三个文件,每个文件里面都是对应RNA的gene name。
- 使用grep命令遍历文件中的gene name,在gtf文件中查找相应的内容。
注意:这里需要使用参数-w来进行精确匹配,不然会得到许多包含该字符串的内容
(说多了都是泪啊/(ㄒoㄒ)/~~)
1.2gtf转化为bed文件
在这里,我们了解一下gtf文件和bed文件有什么特点
NGS基础 - GTF/GFF文件格式解读和转换
言归正传,能够将gtf文件转化为bed的文件有很多,上面那篇推文中说到的Cufflink可以选择,在这里我使用的是bedops中gtf2bed
这个命令。
2. bed提取真实数据bam
samtools view -hb input.bam -L 制作的.bed > output.bam
经过此操作便可以得到真实样本中不同RNA种类所对应的bam文件。
因为这里我们需要统计的是clean reads的长度,并不是fragment的长度,因此有些现成的工具我们还不能使用。
3. bam转化fastq
根据我们生成的bam在clean reads的fastq文件中将不同类型RNA对应的fq提取出来。
samtools view input.bam | awk '{print $1}' | sort | uniq > output_name.list
获取相应的reads序号。
seqtk subseq clean_reads.1.fq.gz output_name.list > output.fq
这里我在提取的时候报错,说什么格式不对,我就在每一行后面添了/1。
sed -i 's/$/\/1/g' *list
之后就可以提取序列,得到不同RNA对应的fq文件。
4. 统计fq文件中reads长度
有工具可以统计fq文件中reads的长度(忘记了,命令里有个fax…啥的,之后想起来再补上)。
但我不理解为什么我当时统计出来什么都没有,没办法自己写咯。
for rnatype in lncRNA miRNA miscRNA Mt_tRNA Mt_rRNA other protein_coding pseudo rRNA scRNA snoRNA snRNA tRNA
do
awk 'NR>1 && (NR-2)%4==0{print length}' $rnatype.fq >> ${rnatype}_length.txt
done