该页面作为开题答辩的补充页，对答辩研究进度部分信息进行补充，只展示部分关键数据和代码，更细节部分请问我
下划线文字附带超链接，请自行浏览。只展示关键分析内容和结果，最新进展请以本人的嘴为准。
1.表型数据分析
1.1统计性直方图
根系性状符合正态分布，说明其为受微效基因调控的数量性状。
1.2 统计性数据分析
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-011.3相关性分析
两种呈现方式
2.基因型数据预处理
2.1基因型数据质控
数据质控是为了
- 增加计算结果的准群性
- 提高运行计算的速度和效率
2.1.1基因型数据缺失质控
在git环境下基于PLINK程序进行，共818169 个SNP位点 and 476 个个体 pass filters and QC.
plink --bfile 2 --geno 0.1 --make-bed --mind 0.15 --out 3_missing通过R语言可视化SNP位点缺失率和个体缺失率
2.1.2maf质控
统计次等基因频率MAF，剔除次等基因频率过低snp位点以提高结果正确性和运算速度
在git环境下基于PLINK程序进行（如上）
通过R语言可视化SNP位点MAF分布
2.1.3杂合率质控
在git环境下基于PLINK程序进行,818169 variants and 476 people pass filters and QC.
plink --bfile 3 --het --out het_check2.1.4 数据导出
剩余的质控对于我们所用的植物SNP数据并无太多影响，基因型数据质控到此为止，通过plink封装数据后导出
2.2 群体结构的验证
这一部分我们的计算与已有的研究虽有差距，但是大体上相同，是一次尝试。
2.2.1祖先成分的分析
我们的组先成分分析基于部署在腾讯云的云服务器进行，基于的经典的祖先成分分析程序Admixture进行
我们的分析过程如下
- 基于数据降维和admixture对于各亚群数量（K的值）计算得到的CV值，确定K=3-4，这与已有的研究相同
- 基于得到的K值进行祖先成分的计算
- 通过R语言可视化分析结果
我们的祖先成分图没经过精细排序，看起来有点乱，但是可以看出被划分为了3个主要祖先成分，与已有的研究相符。2.2.2群体结构的解析
群体结构的解析我们基于tassel 5.0构建邻接树，导出数据后通过R语言绘制环形邻接树，通过FigTree.v1.4.4对不同分支进行上色和美化
主要的分支状况与该群体的研究类似，至此我们确定了已有研究的群体结构数据还是很可行的（hhh）2.2.3 LD距离的计算
这也是一次尝试，与已有的研究差距不是很大，所以我们的研究选择50 kb寻找候选基因没什么问题。
ld <- read.delim(file.choose(),stringsAsFactors = FALSE,header=TRUE, sep = "\t")
str(ld)
##remove sites that have NaN for distance or r2
ld_sub <- ld[ld$R.2 != "NaN",]
ld_sub$dist <- as.numeric(ld_sub$Dist_bp)
ld_sub2 <- ld_sub[ld_sub$dist != "NaN",]
ld_sub2$rsq <- ld_sub2$R.2
file <- ld_sub2[,c(1,2,7,8,15:19)]
# C values range from about 0.5 to 2, start with 0.1
Cstart <- c(C=0.1)
# fit a non linear model using the arbitrary C value,
# N is the number of the genotypes that have the SNP site
modelC <- nls(rsq ~ ( (10+C*dist)/( (2+C*dist) * (11+C*dist) ) ) *
( 1+( (3+C*dist) * (12+12*C*dist+(C*dist)^2) ) / ( 2*N*(2+C*dist) * (11+C*dist) ) ),
data=file, start=Cstart, control=nls.control(maxiter=100))
# extract rho, the recombination parameter in 4Nr
rho <- summary(modelC)$parameters[1]
# feed in the new value of rho to obtain LD values adjusted for their distances along the chromosome/genome
newrsq <- ( (10+rho*file$dist) / ( (2+rho*file$dist) * (11+rho*file$dist) ) ) *
( 1 + ( (3+rho * file$dist) * (12+12*rho*file$dist + (rho*file$dist)^2) ) /
(2*file$N*(2+rho*file$dist) * (11+rho*file$dist) ) )
newfile <- data.frame(file$dist, newrsq)
maxld <- max(newfile$newrsq,na.rm=TRUE) #using max LD value from adjusted data
halfdecay = maxld*0.5
halfdecaydist <- newfile$file.dist[which.min(abs(newfile$newrsq-halfdecay))]
newfile <- newfile[order(newfile$file.dist),]
# plotting the values
tiff("LD3.tif",width = 800,height = 800)
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 4, 0, 0))
plot(file$dist, file$rsq, pch=".", cex=8, xlab="Distance (bp)", ylab=expression(LD ~ (r^2)), col="grey")
lines(newfile$file.dist, newfile$newrsq, col="red", lwd=2)
abline(h=0.2, col="blue") # if you need to add horizental line
abline(v=halfdecaydist, col="green")
mtext(round(halfdecaydist,2), side=1, line=0.05, at=halfdecaydist, cex=0.75, col="green")
dev.off()
一般我们是选取0.2为阈值，图上标注的是降到一半0.5的数值为5 kb，降低到0.2约为50kb。我们之后找候选基因就用50kb了。
3.1PCA分析
3.1.1基于样本的PCA
通过R语言进行分析
碎石图，表示各个PCA维度的解释率，两个维度已经不错了，我们就选择2维度PCA
PCA分析，两个性状从根系8个性状中分离出来
3.1.2 基于性状的PCA
AREA与RTP性状分离出来，AREA解释大部分的变异
4.1亚群数据分析
亚群数据基于我们所得到的群体结构数据，基于此对421个自交系进行亚群分类，亚群间的多重比较基于SPSS进行，结果通过R语言“boxplot”进行可视化
5.1GWAS分析
5.1.1GWAS分析过程
GWAS的分析基于tassel 5.0进行，通过结合群体结构数据、kinship数据和基因型数据计算。选定模型为混合线性模型MLM。
阈值的确定基于广泛接受的方法确定，即以1除以通过Genetic Type I error方法计算得到独立的有效SNP个数（490547）确定阈值p为2.039E-06，以此判定SNP标记与目标性状是否关。
为了数据的美观我们通过R语言绘制了环形曼哈顿图
5.1.2候选基因的选择
候选基因的选择我们基于NCBI与MaizeGDB进行，选择SNP位点前后50 kb范围筛选候选基因。同时基于TARI我们寻找了该基因的拟南芥同源基因并进行注释，以对无基因注释的基因功能进行一定地暗示。
为了进行大批量SNP数据的注释，我们通过R语言进行基因映射和注释,其原理为先下载基因轨道数据和注释数据后根据选取100kb区段的物理位置直接进行映射
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 = "、")
}我们共寻找到好多（189）个候选基因。
6.后续分析
候选我们将通过共表达网络进一步筛选这189个候选基因，并对其中高可行度的基因进行连锁分析、单倍体分析、基因表达分析等。