找到了一份種群基因組學(xué)數(shù)據(jù)分析的教程，原文用的數(shù)據(jù)是2015年發(fā)表在science上的一篇論文Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake。這份教程利用這篇文章的數(shù)據(jù)分析了部分內(nèi)容。

教程地址

https://github.com/singhal/popgen_tutorial/blob/master/pop_gen_tutorial.rst

vcf文件下載鏈接

http://datadryad.org/bitstream/handle/10255/dryad.101389/Massoko_Dryad_VCF_final.vcf.gz

爭取把教程中的內(nèi)容都重復(fù)完，然后多看幾遍原論文（高水平論文看起來還真有些吃力呢?。?/p>

下載數(shù)據(jù)

wgethttp://datadryad.org/bitstream/handle/10255/dryad.101389/Massoko_Dryad_VCF_final.vcf.gz

為了減小計(jì)算壓力，教程中的處理方式是只保留36個樣本（正常數(shù)據(jù)中好像是有146個樣本，解壓出來的vcf文件有11G），并且刪除了inde，只保留snp位點(diǎn)。但是原文中保留的36個個體的文本文件inds_to_keep.txt我現(xiàn)在找不到，需要自己重新構(gòu)造一份需要保留的個體的樣本名。處理方式是：

首先使用bcftools工具將所有的樣本名重定向到一個文件里

bcftoolsquery-lMassoko_Dryad_VCF_final.vcf.gz>inds_to_keep.txt

我選擇的是每個群體保留六個樣本（樣本前綴名一直我就認(rèn)為他們是來自同一個群體），最后我保留了38個個體

這一步大家可以自行進(jìn)行處理或者給我留言獲得inds_to_keep.txt文件。

提取指定的樣本并刪除indel

vcftools--gzvcf?Massoko_Dryad_VCF_final.vcf.gz?--keep?inds_to_keep.txt?--stdout?--recode?--recode-INFO-all--remove-indels?|?bgzip??>?Massoko_Dryad_VCF_final_subset_noIndels.vcf.gz

為了減小計(jì)算壓力，進(jìn)一步對文件進(jìn)行處理（這一步使用到的兩個參數(shù)自己還不太明白是什么意思，這一步完全照搬原教程）

vcftools--gzvcf?Massoko_Dryad_VCF_final_subset_noIndels.vcf.gz?--maf?0.05?--max-maf?0.95?--stdout?--recode?--recode-INFO-all|?bgzip?>?Massoko_Dryad_VCF_final_subset_noIndels_maf05.vcf.gz

vcftools--gzvcf?Massoko_Dryad_VCF_final_subset_noIndels_maf05.vcf.gz?--thin?1000?--stdout?--recode?--recode-INFO-all|?bgzip?>?Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf.gz

這里不明白的參數(shù)

--maf

--max-maf通常會設(shè)置最小等位基因頻率來過濾vcf文件，但這里設(shè)置最大等位基因頻率是什么意思？

--thin 1000

接下來計(jì)算兩個不同群體的核苷酸多樣性

獲得兩個不同群體所有的樣本名，存入文件中

bcftools?query?-l?Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf?|grep"littoral">?littoral.txt

bcftools?query?-l?Massoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf?|grep"benthic">?benthic.txt

計(jì)算群體核苷酸多樣性

vcftools--vcfMassoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf--keeplittoral.txt--window-pi100000--outlittoral_pi

vcftools--vcfMassoko_Dryad_VCF_final_subset_noIndels_maf05_thinned1K.vcf--keepbenthic.txt--window-pi100000--outbenthic_pi

--window-pi?指定窗口的長度

--out?指定輸出文件的前綴名

將結(jié)果文件導(dǎo)出，使用ggplot2做折線圖和箱線圖

箱線圖

bb<-read.table("../../vcf_handling/Fish_Populations/benthic_pi.windowed.pi",header=T)

ll<-read.table("../../vcf_handling/Fish_Populations/littoral_pi.windowed.pi",header=T)

dim(bb)

head(bb)

bb$indiv<-"benthic"

dim(ll)

head(ll)

ll$indiv<-"littoral"

df<-rbind(bb,ll)

colnames(df)

dim(df)

library(ggplot2)

options(scipen=200)

ggplot(df,aes(x=indiv,y=PI,fill=indiv))+

geom_boxplot()+theme_bw()

image.png

折線圖

ggplot(df,aes(x=BIN_START,y=PI,group=indiv,color=indiv))+

geom_line()+theme_bw()+

theme(legend.position?="top",

legend.title?=?element_blank(),

axis.text.y?=?element_blank(),

axis.ticks.y?=?element_blank(),

axis.title?=?element_blank())

image.png