该页面作为开题答辩的补充页，对答辩部分信息进行补充，只展示部分关键数据和代码，更细节部分请问我（手下留情啊~）
下划线文字附带超链接，请自行浏览。
1.1试验材料
本实验所使用的基因型数据来自于MaizeGO，数据类型为SNP数据。GWAS所使用的群体结构来自于该自交群体已有的研究
（Yang, X., Gao, S., Xu, S. et al. Characterization of a global germplasm collection and its potential utilization for analysis of complex quantitative traits in maize. Mol Breeding 28, 511–526 (2011)).Spring。
群体结构文件请从http://www.maizego.org/download/M.pdf 下载，该群体共527个自交系，我们的实验基于其中的421个自交系进行。
基因型数据结构为SNP数据，群体结构数据（图1）基于对于玉米群体的SS、NSS、TST亚群划分进行
1.2根系表型的测定
该部分实验基于导师2015年海南已有的实验。具体方式为于成熟期选取田间长势均匀一致的3-5株进行玉米根系的挖掘，挖掘深度为30 cm。先用抖土法抖掉根系附着在根系上的大部分土壤，然后将根系浸泡在混有洗涤剂的水池中，最后再用压力可调的高压清洗机（亿力牌4630C-120C，中国）清洗掉附着在根系上的残留土壤。待根系清洗干净后，将根系转移至光源稳定的摄影棚（DEEP牌80 cm型，中国）中对根系的二维图像的进行采集。采集根系图像的数码相机为SONY-5100L型微单相机（日本）。利用软件DIRT对二维根系图像进行高通量定量化分析。
DIRT的详细使用方法，请点击本站点的另一篇文章——植物根系性状的测定。
1.3表型数据的分析
表型数据分析基于R 4.03进行。
（1）统计性分析
统计性分析基于R语言，对数据进行直方图绘制，并拟合正态曲线
a<-read.table("2015-phenotypic date.txt",sep = "\t",header = T)
#加载包
library(ggplot2)
library(patchwork)
#绘图
#area
a<-a[,c(1,6:8,5,4,2,3,9)]
par(mfrow = c(2, 4))
for (i in 2:9) {
f<-a[,i]
hist(f,prob=T,col="grey",xlab = "",main="")
xfit<-seq(min(f),max(f),length=40)
yfit<-dnorm(xfit,mean(f),sd(f))
lines(xfit,yfit,col="red",lwd=3)
}
(2)统计性分析
统计性分析基于R语言进行，对于其常规统计性结果如均值、CV进行计算
a<-read.table("2015-phenotypic date.txt",sep = "\t",header = T)
library(pastecs)
stat.desc(a[,-1])
c<-stat.desc(a[,-1])
> c
AREA AVG_DENSITY WIDTH_MED WIDTH_MAX ANG_TOP ANG_BTM SKL_WIDTH RTP_COUNT
nbr.val 4.210000e+02 421.0000000 4.210000e+02 4.210000e+02 4.210000e+02 4.210000e+02 4.210000e+02 4.210000e+02
nbr.null 0.000000e+00 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
nbr.na 0.000000e+00 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
min 4.828880e+03 1.4787060 6.819370e+01 7.925440e+01 4.910777e+01 2.978805e+01 8.961045e+01 1.266667e+02
max 2.213200e+04 29.2980000 1.905613e+02 2.427560e+02 7.172780e+01 5.941517e+01 2.529947e+02 1.044000e+03
range 1.730312e+04 27.8192940 1.223676e+02 1.635016e+02 2.262003e+01 2.962712e+01 1.633842e+02 9.173333e+02
sum 5.144396e+06 3121.2662722 4.849089e+04 6.266166e+04 2.625826e+04 1.891334e+04 6.532830e+04 1.864558e+05
median 1.196543e+04 6.3614200 1.132773e+02 1.456183e+02 6.265830e+01 4.482257e+01 1.534190e+02 4.296667e+02
mean 1.221947e+04 7.4139341 1.151802e+02 1.488400e+02 6.237116e+01 4.492481e+01 1.551741e+02 4.428878e+02
SE.mean 1.552185e+02 0.2023009 1.056573e+00 1.420116e+00 2.080944e-01 2.751109e-01 1.446403e+00 7.221837e+00
CI.mean.0.95 3.051018e+02 0.3976483 2.076830e+00 2.791419e+00 4.090363e-01 5.407658e-01 2.843091e+00 1.419545e+01
var 1.014306e+07 17.2296944 4.699821e+02 8.490426e+02 1.823069e+01 3.186381e+01 8.807668e+02 2.195722e+04
std.dev 3.184817e+03 4.1508667 2.167907e+01 2.913834e+01 4.269741e+00 5.644804e+00 2.967772e+01 1.481797e+02
coef.var 2.606346e-01 0.5598737 1.882186e-01 1.957695e-01 6.845697e-02 1.256500e-01 1.912543e-01 3.345761e-01(3)相关性分析
相关性分析基于R语言，通过皮尔逊相关计算其相关性矩阵和P值矩阵，通过corrplot展示相关性矩阵。有下图两种展示方式。
library(Hmisc)
library(PerformanceAnalytics)
a<-read.table("2015-phenotypic date.txt",sep = "\t",header = T)
#相关性分析
d<-a[,-1]
d<-as.matrix(d)
res<-rcorr(d)
#绘制相关性分析图
png(filename = "ralation.png")
chart.Correlation(d,histogram = T,method = c("pearson"))
dev.off
#绘制相关性热图
library(corrplot)
data<-cor(d)
pdf("ralation cor.pdf")
corrplot(data,type = "upper",order="hclust",tl.col = "black",tl.srt = 45,tl.cex = .8)
dev.off
1.4 GWAS分析
基于Tassel 5.0软件进行。通过SNP数据和群体结构以及计算得到的kinship矩阵进行基于混合线性模型（MLM）的全基因组关联分析。关于Tassel 5的详细操作方法，请见本站点的另一篇文章——通过Tassel进行GWAS分析。
在此分析之前，我们尝试通过独立计算群体结构以验证群体结构数据的准确性。我们进行了如下尝试：
- 先进行数据降维和PCA分析得到该群体适宜的亚群数目，K=3-4，与已有研究相符合。基于云服务器，通过ADMIXTURE: fast ancestry estimation进行祖先成分的分析。这结果已有的研究大致相符合。
- 通过Tassel 5分析群体结构，分析结构绘制邻接树与已有的研究进行对比。与已有的研究大致相符合。
因此我们使用已有研究的群体结构数据进行GWAS分析。library(stringr) library(ggtree) library(treeio) pdf("tree222.pdf") tree<-read.tree("tree.nwk") ggtree(tree,branch.length = 'none',layout = 'circular')+ geom_tiplab(size=1) dev.off
GWAS得到的结果我们通过R语言进行绘制曼哈顿图和qq图
#Manhattan plot
setwd("H:/GWAS-MLM/GWAS data")
#
library(qqman)
#ANG_BTM
png("ANG_BTM.png",width = 1107,height = 556)
a<-read.table("ANG_BTM.txt",sep = '\t',header = T)
a<-a[-1,]
manhattan(a,chr="Chr",bp="Pos",p="p",snp = "Marker",
main="ANG_BTM",
suggestiveline = 5.691,
genomewideline = FALSE
)
dev.off
#剩余7个性状绘制类似library(qqman)
lay1<-lay_new(mat = matrix(1:8, ncol = 4))
lay_set(lay1)
#ANG_BTM
a<-read.table("ANG_BTM.txt",sep = '\t',header = T)
a<-a[-1,]
qq(a$p,main="ANG_BTM")
#剩余7个性状绘制类似
为了美观，我们绘制了多组性状的环状曼哈顿图
data<-read.table("datanew.txt",sep = '\t',header = T)
library(CMplot)
a<-c(1:3,5,9,10,6,11,8,7,4)
datanew<-data[a]
CMplot(datanew,type="p",plot.type="c",r=0.4,col=c("grey30","grey60"),chr.labels=paste("Chr",c(1:10),sep=""),
threshold=c(2.039e-6),cir.chr.h=1.5,amplify=TRUE,
threshold.lty=c(1,2),threshold.col=c("blue"),
signal.line=1,signal.col=c("red","green"),
chr.den.col=NULL,
ylim = c(0,8),
signal.cex = 0.6,
cir.legend = T,
cir.band = 1.5,
lwd.axis = 0,
bin.size=1e4,outward=FALSE,file="jpg",
memo="",dpi=300,file.output=TRUE,verbose=TRUE,width=10,height=10)下图只展示部分图片
我们在显著SNP位点的前后50kb范围寻找候选基因，玉米基因组的轨道文件源于Maize GDB
我们通过R语言进行批量的候选基因寻找
library(readxl)
snp<-read_xlsx("地下部LOD大于5.xlsx",sheet = 1)
gene<-read.table("ZmB73_5b_FGS-v2_genepos.txt",header = F,sep = "\t")
#重组gene数据
names(gene)<-c("Marker","chr","start","end")
#重组SNP数据
snp<-data.frame(snp[,c(1:3,6,4)])
#确定50kb范围
snp$pos.l<-snp$Pos-50000
snp$pos.r<-snp$Pos+50000
#开始寻找基因，内外双循环结构
for (i in 1:nrow(snp)) {
rowname<-c() #存储基因行名的暂存向量
for (f in 1:nrow(gene)) {
if(snp[i,3]==gene[f,2]|gene[f,3]>=snp[i,6]|gene[f,4]<=snp[i,7]){
rowname<-c(rowname,f)
}
}
snp$allgene<-paste(gene[rowname,1],collapse = "、")
}我们共找到显著SNP关联的189个候选基因
1.5共表达网络的构建和分析
我们通过GWAS得到的189个候选基因作为共表达网络的输入基因列表，充当差异表达基因的角色。我们计划基于Camoco对我们GWAS筛选的显著SNP与前后50kb范围内的基因进行映射。计算基因间的共表达，并标记其显著性。计算出的结果将通过Cytoscape软件展示结果。
我们进行共表达网络的构建有如下目的
- 剔除在大部分样品中不表达的候选基因，进一步筛选候选基因。
- 根据共表达网络的共表达状况，我们挑选优先度较高的基因进行进一步的分析。
1.6基因表达量的测定
这一部分实验基于导师2019年已有的实验，通过RNA-seq进行。具体步骤为通过Trizol试剂盒进行总RNA的提取后进行评估RNA完整性，并通过Illumina HiSeq进行测序。测序得到的数据通过HISAT2进行序列对比，最后通过StringTie将基于表达水平组装为FPKM值，即将测序仪得到的reads数据转化为每千个碱基的转录每百万映射读取。我认为这是一类标准化过程，以比对基因表达的差异。
