【毕业开题】开题试验设计信息详情

该页面作为开题答辩的补充页，对答辩部分信息进行补充，只展示部分关键数据和代码，更细节部分请问我（手下留情啊~）

下划线文字附带超链接，请自行浏览。

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数据转化为每千个碱基的转录每百万映射读取。我认为这是一类标准化过程，以比对基因表达的差异。