简化基因组群体分析

1.stack寻找位点
1)数据预处理
数据下机经过拆分去接头获得原始数据

图片.png

由于拆分数据用的引物街头长度不一，获得的原始数据的长度也有差异，这对后续stack分析有影响，这里为了减少后续分析的麻烦截取reads的前135bp
使用工具seqtk

seqtk trimfq -L 135 ../data/FNA-7.2.fq.gz > FNA-7.2.fq

2）stack运行流程
由于没有近缘的参考基因组，采用无参流程。该流程参考stack官网步骤（http://catchenlab.life.illinois.edu/stacks/manual/#popmap）。

#!/bin/bash

src=$HOME/research/project

files=”sample_01
sample_02
sample_03”

#
# Build loci de novo in each sample for the single-end reads only. If paired-end reads are available, 
# they will be integrated in a later stage (tsv2bam stage).
# This loop will run ustacks on each sample, e.g.
#   ustacks -f ./samples/sample_01.1.fq.gz -o ./stacks -i 1 --name sample_01 -M 4 -p 8
#
id=1
for sample in $files
do
    ustacks -f $src/samples/${sample}.1.fq.gz -o $src/stacks -i $id --name $sample -M 4 -p 8
    let "id+=1"
done

# 
# Build the catalog of loci available in the metapopulation from the samples contained
# in the population map. To build the catalog from a subset of individuals, supply
# a separate population map only containing those samples.
#
cstacks -n 6 -P $src/stacks/ -M $src/popmaps/popmap -p 8

#
# Run sstacks. Match all samples supplied in the population map against the catalog.
#
sstacks -P $src/stacks/ -M $src/popmaps/popmap -p 8

#
# Run tsv2bam to transpose the data so it is stored by locus, instead of by sample. We will include
# paired-end reads using tsv2bam. tsv2bam expects the paired read files to be in the samples
# directory and they should be named consistently with the single-end reads,
# e.g. sample_01.1.fq.gz and sample_01.2.fq.gz, which is how process_radtags will output them.
#
tsv2bam -P $src/stacks/ -M $src/popmaps/popmap --pe-reads-dir $src/samples -t 8

#
# Run gstacks: build a paired-end contig from the metapopulation data (if paired-reads provided),
# align reads per sample, call variant sites in the population, genotypes in each individual.
#
gstacks -P $src/stacks/ -M $src/popmaps/popmap -t 8

#
# Run populations. Calculate Hardy-Weinberg deviation, population statistics, f-statistics
# export several output files.
#
#populations -P $src/stacks/ -M $src/popmaps/popmap -r 0.65 --vcf --genepop --structure --fstats --hwe -t 8
populations -P run1/ -M popmap_sample.tsv -r 0.8 --vcf --genepop --structure --hwe -t 30 --fasta-samples --plink --phylip_var_all  -O test #最后一步按照自己的需求获取相应的文件

1-1.位点过滤
得到原始的vcf文件后，一般位点会比较多，可能会包含许多假阳性位点等。所以进行一步过滤操作。

vcftools --gzvcf populations.snps.vcf.gz --max-missing 0.5 --mac 3  --recode --recode-INFO-all --out populations_filt 
###最大缺失率50%，等位基因次要基因型频数最少3

获得过滤后的vcf文件进行后续分析。

2.PCA分析
参考（https://zhuanlan.zhihu.com/p/45950835）

library(gdsfmt)
library(SNPRelate)
vcf.t <- "study.vcf"
snpgdsVCF2GDS(vcf.t, "new_geno.gds", method = "biallelic.only")
genofile <- snpgdsOpen("new_geno.gds")
pca <- snpgdsPCA(genofile,autosome.only=FALSE)
pc.percent <- pca$varprop * 100
tab <- data.frame(sample.id = pca$sample.id, EV1 = pca$eigenvect[,1], EV2=pca$eigenvect[,2],EV3=pca$eigenvect[,3],EV4=pca$eigenvect[,4],stringsAsFactors =F)
plot(tab$EV2, tab$EV1, xlab="PC2", ylab="PC1")

图片.png

3.进化树分析

set.seed(100)
ibs.hc <- snpgdsHCluster(snpgdsIBS(genofile, num.thread=2,autosome.only = FALSE))
rv <- snpgdsCutTree(ibs.hc)
plot(rv$dendrogram, leaflab="perpendicular", main="Tree")

这里有个地方很奇怪，snpgdsIBS计算的时候矩阵会出现NaN，导致后面snpgdsHCluster聚类不成功，现在不知道什么原因，只能把NaN替换成0

b<-snpgdsIBS(genofile, num.thread=2,autosome.only = FALSE)
b$ibs[is.na(b$ibs)]<-0
ibs.hc <-snpgdsHCluster(b)

图片.png
4.Structure分析
vcf转bed参考（https://www.jianshu.com/p/dc82fcbe3cda）第一列要修改一下
使用admixture软件,输入plink格式的bed文件。
跑之前要对数据先处理一下

plink --vcf populations_filt.recode.vcf --recode --out Arab  --allow-extra-chr
cat  populations.plink.map | sed '1d' | awk '{print "1\t"$2"\t"$3"\t"$4}' > populations.plink1.map #第一列un软件识别不了，改成1
plink --make-bed --ped  populations.plink.ped --map populations.plink1.map

admixture plink.bed 2 
admixture plink.bed 4 

最终获得两个文件plink.2.P，plink.2.Q
画图

tab<-read.table('plink.2.Q')
barplot(t(as.matrix(tab)),col=rainbow(2),xlab="Individual #",ylab="Ancestry",border=NA)

图片.png