拆分bam统计clean reads长度

因为需要统计不同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文件

1. 获取RNA type和gene name的对应关系，这里是根据ensembl网站上下载的数据先得到gene name对应的RNA类型。这里因为工作需要大概分成了十三种，可以得十三个文件，每个文件里面都是对应RNA的gene name。
2. 使用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