bam文件的可视化（测序深度） | IGV

写在前面

没有系统地学过生信，当时靠着导师的“威逼”，在学长的指导下学了perl，然后开始了边学边用的生信路程。

回头去看，才意识到基础有多不扎实。所以现在但凡遇到+解决一点小问题，都尽量归纳整理出来。

之前用IGV去看比对情况，直接导入bam文件(一个外显子的比对文件6G左右），结果笔记本就卡顿了，特别不顺。当时因为不着急，也就没想过要去了解变通方案。

最近因为分析meRIP-seq的数据，意外在过程中发现，可以将bam转为tdf文件。tdf文件很小，再也不用担心IGV卡顿了。但是tdf文件只能反映基因组每个区域的测序深度，无法看到具体的比对情况，适合用来check找到的peak或者CNV。

实战

#安装igvtools
conda install igvtools

#remove duplicated reads
samtools rmdup -s sample.bam sample.rmdup.bam 
samtools index sample.rmdup.bam

#自定义genome的chrom.sizes文件 （根据bam的header获得），这里需要根据具体的header来去除相应的行，仅保留染色体名称和长度信息
samtools view -H sample.rmdup.bam | awk -F '[\t:]' '{print $3"\t"$5}' | grep -v "coordinate" | grep -v "bwa" | grep -v "SAM" >mm10.chrom.sizes

cp mm10.chrom.sizes /root/miniconda2/share/igvtools-2.3.93-0/genomes/

#从bam生成tdf
igvtools count -z 5 -w 10 -e 0 sample.rmdup.bam sample.tdf mm10 
#-w The window size over which coverage is averaged. Defaults to 25 bp.
#The count command computes average feature density over a specified window size across the genome.
#The input file must be sorted by start position. See the sort command below.
#会调用 /root/miniconda2/share/igvtools-2.3.93-0/genomes/mm10.chrom.sizes 文件，需要chrom名字对应

将tdf文件导入本地版IGV，选择相应的基因组（这里是mm10）,就可以查看了。

Note : 如果发现一片空白，好像啥都没有，不要心慌。可能只是显示了全基因组，信息太多没显示出来，选择其中一条染色体，也就是缩小显示范围，就能看到希望的柱子了。O(∩_∩)O哈哈~