沒有系統(tǒng)地學過生信,當時靠著導師的“威逼”,在學長的指導下學了perl,然后開始了邊學邊用的生信路程。
回頭去看,才意識到基礎有多不扎實。所以現(xiàn)在但凡遇到+解決一點小問題,都盡量歸納整理出來。
之前用IGV去看比對情況,直接導入bam文件(一個外顯子的比對文件6G左右),結(jié)果筆記本就卡頓了,特別不順。當時因為不著急,也就沒想過要去了解變通方案。
可以將bam轉(zhuǎn)為tdf文件。tdf文件很小,再也不用擔心IGV卡頓了。但是tdf文件只能反映基因組每個區(qū)域的測序深度,無法看到具體的比對情況,適合用來check找到的peak或者CNV。
#安裝igvtoolsconda install igvtools#remove duplicated readssamtools rmdup -s sample.bam sample.rmdup.bam samtools index sample.rmdup.bam#自定義genome的chrom.sizes文件 (根據(jù)bam的header獲得),這里需要根據(jù)具體的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.#會調(diào)用 /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 : 如果發(fā)現(xiàn)一片空白,好像啥都沒有,不要心慌。可能只是顯示了全基因組,信息太多沒顯示出來,選擇其中一條染色體,也就是縮小顯示范圍,就能看到希望的柱子了。O(∩_∩)O哈哈~