群体进化分析

一、SNP预处理
必备小知识：
PLINK是一种广泛使用的生物信息学工具，它支持多种文件格式来存储和处理遗传数据。PLINK的文件格式主要分为两类：文本格式和二进制格式。

文本格式：map/ped文件

map文件：存储变异位置信息，包括染色体编号、变异标识符、摩尔位置和碱基对坐标。map文件没有列名，通常使用四列数据表示。

ped文件：记录个体的系谱信息和基因型数据，包括家系ID、个体ID、父母ID、性别和表型值。ped文件效率较低，不建议用于下游分析。

二进制格式：bed/bim/fam文件

bed文件：存储基因型信息，是PLINK中的二元等位基因表。与UCSC Genome Browser的BED格式不同，PLINK的bed文件必须与bim和fam文件一起使用。

bim文件：存储每个遗传变异的相关信息，包括染色体编号、变异标识符、摩尔位置、碱基对坐标和等位基因信息。

fam文件：存储样本家系等信息，包括家系编号、个体编号、父母编号、性别和表型值。

1.VCF格式转为plink文本格式，得到ped和map文件

vcftools \                        # 软件程序
--gzvcf Filter.snp.vcf.gz \       # 指定输入文件（要求vcf格式的gz压缩文件）
--chrom-map chrom_map.txt \       # 指定染色体名的映射文件
--maf 0.05 \                      # 指定次等位基因频率的过滤标准
--plink \                         # 转成plink格式
--out ramie_plink1                # 指定输出文件名的前缀
#chrom_map.txt文件格式
Chr1    1
Chr2    2
Chr3    3
Chr4    4
Chr5    5
Chr6    6
Chr7    7
Chr8    8
Chr9    9
Chr10   10

2.过滤出连锁平衡位点，得到prune.in文件

plink \                         # 软件程序
--threads 4 \                   # 设置最大线程数
--memory 5000 \                 # 设置初始工作区大小（单位：MB）
--file ramie_plink1 \           # 指定输入文件名的前缀
--indep-pairwise 50 5 0.5 \     # 设置窗口大小(kb)、步长大小(ct)和r^2阈值
--out ramie_plink2              # 指定输出文件名的前缀
  1. 根据.in文件从总的plink文件中提取出连锁平衡的位点,得到map和ped文件，用于后续的群体分层分析
plink \                            # 软件程序
--threads 4 \                      # 设置最大线程数
--memory 5000 \                    # 设置初始工作区大小（单位：MB）
--file ramie_plink1 \              # 指定输入文件名的前缀（总的plink文件）
--extract ramie_plink2.prune.in \  # 指定第2步输出的包含连锁平衡位点的文件
--recode \                         # 将二进制文件转换成可读的文本文件
--out ramie_plink3                 # 指定输出文件名的前缀
  1. 可以将得到的文本格式文件转为二进制文件
plink --noweb --file ramie_plink3 --make-bed --out SAMP
#输出结果：.bed .bim .fam文件

二、PCA分析
1.gcta64软件进行PCA分析

gcta64 --bfile ../../1.Plink/SAMP --make-grm --autosome --out SAMP_grm
# --bfile 读取的是plink的二进制文件
# --autosome 是利用常染色体上的所有SNP信息, 这个不能省略
# --make-grm生成grm矩阵
# --out 生成前缀名
gcta64 --grm SAMP_grm --pca 3 --out SAMP_pca
# --grm grm文件
# --pca PCA的数目为3
# --out 结果输出文件, .eigenval  和.eigenvec文件

2.R里画图

date <- read.table("SAMP_pca.eigenvec",header=F)
head(date)
names(date) = c("Fid","ID","PC1","PC2","PC3")
data2<-read.table("SAMP_pca.eigenval",header=F)###读取各主成分贡献度数据
data3<-data2[c(1:(nrow(data2)-1)),]
PC1contri<-round(data3[1]*100/sum(data3),digits=2)###PC1的贡献率
PC2contri<-round(data3[2]*100/sum(data3),digits=2)###PC2的贡献率
PC3contri<-round(data3[3]*100/sum(data3),digits=2)###PC3的贡献率
xlab<-paste("PC1(",PC1contri,"%)",sep="")###设置X轴坐标标签
ylab<-paste("PC2(",PC2contri,"%)",sep="")###设置Y轴坐标标签
zlab <- paste("pc3(",PC3contri,"%)",sep="")
library(ggplot2)
library(scatterplot3d)
p <- ggplot(date,aes(PC1,PC2)) + geom_point() +labs(x=xlab,y=ylab) +theme_bw()
pca_3d <- scatterplot3d(date$PC1,date$PC2,date$PC3,color = rep(c("#00AFBB"),each=562),pch = 15,lty.hide = 2,
                        xlab = xlab,ylab = ylab,zlab = zlab)
p
pca_3d

三、进化树分析

  1. 计算距离矩阵：
plink \                                  # 软件程序
--file ramie_plink3  \                   # 指定plink格式输入文件前缀
--distance square 1-ibs flat-missing \   # 计算距离矩阵
--out tree                               # 指定输出文件名前缀
  1. 将距离矩阵文件转换成phylip格式的文件
perl mdist2phylip.pl \         # 执行perl脚本
tree.mdist \                   # 指定第1步输出的距离矩阵文件
tree.mdist.id \                # 指定第1步输出的样本名文件
tree.phy.dist                  # 指定输出的文件名
  1. 使用PHYLIP软件的neighbor程序进行NJ法建树，得到outtree文件
neighbor                     # 打开软件程序（交互式程序）或者直接输入绝对路径如：/data8/home/jclong1/Liluli/Evolution_and_Gwas/software/phylip-3.697/exe/neighbor
tree.phy.dist                # 指定第2步输出的phylip格式的距离矩阵文件
Y                            # 确认以上信息
4.outtree文件在itol上美化

进化树美化操作详解- iTOL - 组学大讲堂问答社区 (omicsclass.com)

四、结构图

  1. admixture软件分析
for i in {1..10};
do admixture --cv filter.snp.bed $i >> log.txt
done
# --cv 输入文件为plink二进制bed文件
#$i  K值
# >>log.txt 从1到10运行的out结果都整合到log.txt中，以便grep提取cv行

2.提取CV_error

grep CV log.txt | awk -F ':' '{print NR"\t"$2}' | sed '1i\K\tCV_error' >> CV_for_plot.txt
##  参数
awk -F ':'  以：分隔
      NR  当前行号
      "\t"   空格
      $2   第二列
sed '1i\K\tCV_error' 表示给第一行添加K 空格 CV_error
sed '5i\内容' 表示给第五行添加内容

3.R中画CV_error确定最佳K值

library（ggplot2）
mydata <- read.table("CV_for_plot.txt",header=T，sep="\t")
qplot(x=K,y=CV_error,color=I（'black'）,data=mydata)+geom_line(color="red",lwd=1)
#CV error最小的为最佳K值
  1. 绘制STRUCTURE图
tbl=read.table("filter.snp.3.Q") 
barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA)

五、连锁不平衡衰减（LDdecay）
1.计算整体的LD并绘图

./PopLDdecay/bin/PopLDdecay -InVCF /data8/home/jclong1/Liluli/Evolution_and_Gwas/Rawdata/SAMP562.final.20241106.vcf.gz -OutStat LDdecay
perl PopLDdecay/bin/Plot_OnePop.pl -inFile LDdecay.stat.gz  --output Fig
  1. 计算每个亚群的LD并绘图
~/PopLDdecay/bin/PopLDdecay -InVCF Filter.snp.vcf -SubPop 1-population.txt -OutStat p1.stat -MaxDist 1000
~/PopLDdecay/bin/PopLDdecay -InVCF Filter.snp.vcf -SubPop 2-population.txt -OutStat p2.stat -MaxDist 1000
~/PopLDdecay/bin/PopLDdecay -InVCF Filter.snp.vcf -SubPop 3-population.txt -OutStat p3.stat -MaxDist 1000
~/PopLDdecay/bin/PopLDdecay -InVCF Filter.snp.vcf -OutStat all.stat -MaxDist 2000
#-SubPop 群体ID信息，一行一个ID
#-OutStat 输出文件
#-MaxDist 两个SNP间的最大距离，默认是300kb,最大不超过5mb
~/PopLDdecay/bin/Plot_MultiPop.pl -inList draw.list -output draw -bin1 5000 -bin2 50000
# -bin1 -bin2设置的大一点，曲线会更平滑!
# draw.list 包含上一步生成的文件的绝对路径 以及去掉后缀的文件名
/home/ug0797/data/test-fst/test2/p1.stat.gz p1
/home/ug0797/data/test-fst/test2/p2.stat.gz p2
/home/ug0797/data/test-fst/test2/p3.stat.gz p3
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成，浏览时请结合常识与多方信息审慎甄别。
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容