bam文件的可视化(测序深度) | IGV

写在前面

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

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

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

可以将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文件:File >>> Load from file
导入基因组:Genomes >>> Load genome from file
导入注释文件:File >>> Load from file。从UCSC上可以很方便的下载到.bed格式的注释文件。
选择相应的基因组(这里是mm10),就可以查看了。

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

    原文作者:fatlady
    原文地址: https://www.jianshu.com/p/c0635eb9ac24
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞