拆分bam统计clean reads长度

因为需要统计不同RNA_type对应的clean reads长度，首先需要获得不同RNA对应的gtf文件，转化为bed文件，再从我们真实数据的bam文件中提取出相应的bam，根据提取的bam从包含有clean reads的fastq文件中提取出clean reads，最后统计clean reads的长度，作图。

#mermaid-svg-UxHpnpIynkCrHpOa {font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:16px;fill:#333;}#mermaid-svg-UxHpnpIynkCrHpOa .error-icon{fill:#552222;}#mermaid-svg-UxHpnpIynkCrHpOa .error-text{fill:#552222;stroke:#552222;}#mermaid-svg-UxHpnpIynkCrHpOa .edge-thickness-normal{stroke-width:2px;}#mermaid-svg-UxHpnpIynkCrHpOa .edge-thickness-thick{stroke-width:3.5px;}#mermaid-svg-UxHpnpIynkCrHpOa .edge-pattern-solid{stroke-dasharray:0;}#mermaid-svg-UxHpnpIynkCrHpOa .edge-pattern-dashed{stroke-dasharray:3;}#mermaid-svg-UxHpnpIynkCrHpOa .edge-pattern-dotted{stroke-dasharray:2;}#mermaid-svg-UxHpnpIynkCrHpOa .marker{fill:#333333;stroke:#333333;}#mermaid-svg-UxHpnpIynkCrHpOa .marker.cross{stroke:#333333;}#mermaid-svg-UxHpnpIynkCrHpOa svg{font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:16px;}#mermaid-svg-UxHpnpIynkCrHpOa .label{font-family:"trebuchet ms",verdana,arial,sans-serif;color:#333;}#mermaid-svg-UxHpnpIynkCrHpOa .cluster-label text{fill:#333;}#mermaid-svg-UxHpnpIynkCrHpOa .cluster-label span{color:#333;}#mermaid-svg-UxHpnpIynkCrHpOa .label text,#mermaid-svg-UxHpnpIynkCrHpOa span{fill:#333;color:#333;}#mermaid-svg-UxHpnpIynkCrHpOa .node rect,#mermaid-svg-UxHpnpIynkCrHpOa .node circle,#mermaid-svg-UxHpnpIynkCrHpOa .node ellipse,#mermaid-svg-UxHpnpIynkCrHpOa .node polygon,#mermaid-svg-UxHpnpIynkCrHpOa .node path{fill:#ECECFF;stroke:#9370DB;stroke-width:1px;}#mermaid-svg-UxHpnpIynkCrHpOa .node .label{text-align:center;}#mermaid-svg-UxHpnpIynkCrHpOa .node.clickable{cursor:pointer;}#mermaid-svg-UxHpnpIynkCrHpOa .arrowheadPath{fill:#333333;}#mermaid-svg-UxHpnpIynkCrHpOa .edgePath .path{stroke:#333333;stroke-width:2.0px;}#mermaid-svg-UxHpnpIynkCrHpOa .flowchart-link{stroke:#333333;fill:none;}#mermaid-svg-UxHpnpIynkCrHpOa .edgeLabel{background-color:#e8e8e8;text-align:center;}#mermaid-svg-UxHpnpIynkCrHpOa .edgeLabel rect{opacity:0.5;background-color:#e8e8e8;fill:#e8e8e8;}#mermaid-svg-UxHpnpIynkCrHpOa .cluster rect{fill:#ffffde;stroke:#aaaa33;stroke-width:1px;}#mermaid-svg-UxHpnpIynkCrHpOa .cluster text{fill:#333;}#mermaid-svg-UxHpnpIynkCrHpOa .cluster span{color:#333;}#mermaid-svg-UxHpnpIynkCrHpOa div.mermaidTooltip{position:absolute;text-align:center;max-width:200px;padding:2px;font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:12px;background:hsl(80, 100%, 96.2745098039%);border:1px solid #aaaa33;border-radius:2px;pointer-events:none;z-index:100;}#mermaid-svg-UxHpnpIynkCrHpOa :root{--mermaid-font-family:"trebuchet ms",verdana,arial,sans-serif;}

gtf

bed

样本bam

样本fq

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