全基因組分析後的vcf突變檔案記錄了四百多萬個位點，前面我們講到了如何把它們註釋到dbSNP資料庫ID，一般來說有註釋的位點也就順便註釋到了基因，所以可以簡單寫一個程式來看看哪些基因的突變位點最多:

cat autochr.highQuali.dbsnp.vcf |perl -alne '{/GENEINFO=(.*?)[;,|]/;$h{$1}++ if $1 }END{print "$_t$h{$_}" foreach keys %h}' >~/tmp/vcf.gene.stat

簡單統計結果如下;

當然，其實並不需要註釋到dbSNP資料後再進行統計，可以直接對vcf檔案進行基因註釋，因為vcf檔案本身就記錄著座標，把vcf檔案按照bed格式稍微轉換一下，就可以用bedtools來進行註釋啦。

首先製備好基因的座標檔案，染色體號，基因定位的起始終止座標即可，比如下面這個SPIN1基因：

可以看到， 有10個突變位點註釋到了這個基因，可以其中只有4個是dbSNP資料庫記錄的，所以最開始統計的基因的突變個數排行不是很準確。

那麼我們就針對這個bedtools closest 結果進行統計吧：

可以看到幾乎每個基因的突變個數都增加了，因為不需要被dbSNP資料庫收錄啦。

再看看基因突變個數的個數的變化：

之前突變個數為1的那些基因有1324個，但是現在只剩下了712個！

如果視覺化在圖上，就能看到條形圖很明顯的右移，不過這不是重點。

FAQ

為什麼有712個基因僅僅發現一個突變呢？是這個基因太保守了嗎？還是這個基因長度太短了？同理，那些突變異常多的基因又有什麼特徵呢？

我選取了那712個只有一個變異位點的基因，還有超過400個變異位點的909個基因。

很明顯，從長度來解釋是一個很好的角度~~

以上的變異位點，都應該改名叫做多型性位點。

文：Jimmy

圖文編輯：吃瓜群眾