札記 - 處理兩個物種多個群體的VCF文件中篩選SNP數據

05-07

最近在幫別人做一個近緣物種基因組間漸滲區域和高分化區域的定位分析。

處理基本的SNP數據時候，居然栽了兩回跟頭。活生生一周內搞定的一個分析，被拖長了一倍的時間。

1，從vcf文件裡面call SNP。

根據vcf文件的specification （https://samtools.github.io/hts-specs/VCFv4.2.pdf）

我主要用到的filter條件有(limited to mono or biallelic sites)：

如果一個位點，只有map到一種鹼基，那麼read depth（DP）>=2，認為是純合子genotype，否則認為是N/N。
如果一個位點，map到兩種鹼基，那麼認定這個位點為雜合子的條件是：minor allele supporting read>=2 (AD)
如果一個位點，map到兩種鹼基，並且minor allele只有1條read 支持，那麼如果major allele read depth>=2，則此位點為純合子，否則認為此位點為N/N。

我本來還用到了GQ和PL欄目的信息，GQ越到，對vcf call出來genotypex信心越大。但是在我查取GQ演算法，想要確定一個cutoff value的時候，發現大家對GQ這個參數褒貶不一。並且有人貼了一些比較奇怪的結果，AD支持ancestral和derived allele的數值都很大，但是GQ小於30。For example （pasted from: https://groups.google.com/forum/#!topic/freebayes/fZIjPPuVEkM）

GT:GQ:DP:RO:QR:AO:QA:GL
1) 0/0: 8.95264 : 24 : 24 : 791 : 0 : 0 : 0,-7.22472,-71.5196

2) 0/0:140.745:44:44:1507:0:0:0,-13.2453,-135.973
3) 0/0:140.745:2:2:70:0:0:0,-0.60206,-6.65
4) 0/0 : 140.745 : 12 : 12 : 408 : 0 : 0 : 0,-3.61236,-37.06
1 only has a GQ of 9 despite having 24 reads all being ref. For comparison, 4 has 12 reads and a GQ of 141.
3 has a GQ of 141 with 2 reads. That cant be right, can it?

The RQ values pretty indicate the MQ and BQ are comparable across samples, so it isnt that.

翻譯一下，第一個snp位點，只map到一種鹼基，有24條read的信息，但是GQ只有8.9。

而第四個位點，有12條read map到一種鹼基，GQ卻高達140.

看google網上論壇下面有人說，他從不用GQ來過濾SNP，而是用DP和AD信息。

於是起初我用了GQ>=30的篩選snp條件，後來刪掉了這個條件，只保留了DP和AD的篩選條件。

2，最少？？個體有數據，才call之為SNP

我的species 1中有6個群體，species 2中有4個群體。

這批filter出來的SNP，後面會參與各個群體自己的多樣性計算，兩兩群體之間的遺傳距離計算等。

期初我設定了這樣的filter條件：要求一個位點，在一個物種中>=2個群體中，滿足>=2個個體有genotype信息。

也就是說，經過初步DP&AD的篩選，得到了每一個個體的genotype數據。那麼一個群體只有包含>=2個個體有genotype信息，才認為是一個有效群體。一個物種只有包含>=2有效群體，這個位點才有效。

這樣篩選完後得到的SNP數據，進入下游的多樣性和分化程度的分析。

我自認為設定了一個較為嚴格的篩選條件，但是在後面計算pairwise群體的genetic distance時候，遇到了一點問題。

比如一個位點，在spp1的pop1中有5個個體的信息，在spp2的pop2中有5個個體信息。那麼在pairwise 群體遺傳距離計算中，這個點對於計算pop1和pop2之間的距離應該有貢獻。

但是由於這個位點不滿足上述filter條件（一個物種只有包含>=2有效群體），那麼這個位點就被扔掉了。

也就是說，我扔掉了太多對於pairwise 群體比較中informative的點。於是我修改了filter SNP條件，只要一個物種中包含>=2個有效個體（無論來自1個or幾個群體），那麼這個點就是有效的SNP位點。

3，後面的這個錯誤實在是太小兒科了

我在對一個位點的不同個體call genotype時候，悲催的寫錯了程序。

主要是

如果一個位點，map到兩種鹼基，那麼認定這個位點為雜合子的條件是：minor allele supporting read>=2 (AD)
如果一個位點，map到兩種鹼基，並且minor allele只有1條read 支持，那麼如果major allele read depth>=2，則此位點為純合子，否則認為此位點為N/N。

這麼兩個條件。我一開始寫程序，只考慮了第一種情況，也就是讀到兩種鹼基，minor allele read count>=2的就認為是雜合子。而對於minor allele read count<2的case，直接當成『N/N輸出了。

後來才意識到，minor allele不合格，接下來是查看major allele的count，看是否滿足call homozygous 的條件。

雖然大大小小數據經過手的也不算少，但是，這次犯的錯誤還是讓我印象深刻。其實前面數據處理有問題，總會在後面分析的某一步發覺不對。比如，我是在挑出高變區，畫樹的時候，發覺樹形奇怪，看回每個個體的fasta文件，發現似乎扔掉了太多的點。進而發現，上述第3個錯誤和第1個錯誤的。

發覺不對，回去看raw data，及時發現問題，解決問題。