写在前面没有系统地学过生信,当时靠着导师的“威逼”,在学长的指导下学了perl,然后开始了边学边用的生信路程。 回头去看,才意识到基础有多不扎实。所以现在但凡遇到+解决一点小问题,都尽量归纳整理出来。 之前用IGV去看比对情况,直接导入bam文件(一个外显子的比对文件6G左右),结果笔记本就卡顿了,特别不顺。当时因为不着急,也就没想过要去了解变通方案。 可以将bam转为tdf文件。tdf文件很小,再也不用担心IGV卡顿了。但是tdf文件只能反映基因组每个区域的测序深度,无法看到具体的比对情况,适合用来check找到的peak或者CNV。 实战#安装igvtoolsconda install igvtools#remove duplicated readssamtools 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.sizescp mm10.chrom.sizes /root/miniconda2/share/igvtools-2.3.93-0/genomes/#从bam生成tdfigvtools 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 Note : 如果发现一片空白,好像啥都没有,不要心慌。可能只是显示了全基因组,信息太多没显示出来,选择其中一条染色体,也就是缩小显示范围,就能看到希望的柱子了。O(∩_∩)O哈哈~ |
|