在分析玉米基因组上的一段区域是否受到选择时，可比较玉米基因组和大刍草基因组在这段区域内的多态性，计算多态性指数pi，可以比较这段区域或者基因区域是否受到选择。这里要用玉米材料的群体基因组变异数据和大刍草群体基因组变异数据，并且大刍草的样本和玉米的样本都是回帖在同一套参考基因组上。去年发表的一篇文章，整理了很多前人已经测序的玉米和大刍草样本，并且根据B73参考基因组版本V5重新比对，生成了1500多份包含不同类型大刍草以及玉米的变异数据。用这套数据查看一些基因区域是否选择，非常有代表性。文章发表在the plant journal 上，题目：A common resequencing-based genetic marker data set for global maize diversity 数据下载来源：

https://ars-usda.app.box.com/v/maizegdb-public/folder/189779501832， 该数据把每个染色体分开存放，使用者可以根据自己想要查看的区域，选择合适的染色体进行下载，数据时VCF格式，操作非常方便。同时要下载作者对这些材料进行详细分类的一个表格，这个表格对每个材料的信息进行了详细的分类，比如是玉米还是大刍草，是哪一种大刍草。

1. 准备两个文本文件

一个是玉米的材料编号，一个材料一行，没有表头，另一个是大刍草的编号，也是一个材料一行，没有表头。然后用vcftool进行计算，可以设置窗口大小，按照窗口计算。这里我下载了第四号染色体的文件， 我的目标是比较第四号染色体上，25300000 到 25400000 区间内是否受到选择。每个人可以根据自己实际的情况，选择区域，如果查看一个基因是否受到选择，通常选择这个基因启动子区域，基因区域以及3‘UTR区域，大概几kb。

#maize.csv 文件部分行

11430

1538

178

18_599

196

200B

2369

29MIBZ2

2FACC

2FADB

2MA22

2MCDB

# dachucao.csv 文件部分行

5A1

5A2

5A3

5A4

5A5

5A6

5A7

5B7

5B8

5B9

5D9

5E10

5E11

5E12

5E7

2. Vcftools 进行计算

vcftools --vcf chr_4_imputed.vcf  --chr chr4  --from-bp 25300000 --to-bp 25400000 --keep dachucao.csv --window-pi 1000 --out dachucao_1000.pi

vcftools --vcf chr_4_imputed.vcf  --chr chr4  --from-bp 25300000 --to-bp 25400000 --keep maize.csv --window-pi 1000 --out maize_1000.pi

zemaz = read.table(file='~/maize_hapmap/maize_1000.pi',header = T)

dachucao = read.table(file='~/maize_hapmap/dachucao_1000.pi',header = T)

zemaz$species ='maize'

dachucao$species ='dachucao'

df_plot = rbind(zemaz,dachucao)

library(ggplot2)

ggplot(df_plot,aes(BIN_START,PI,col=species))+geom_line()

write.table(df_plot,file = '~/maize_hapmap/pi_window_1000.txt',row.names = F,quote = F,sep = '\t')