【分享】科研保研小礼包

4.30考研保研分享会补充材料，具体内容如下

常用软件和网站*1
资环PPT模板*4
推荐信模板*1
个人简历模板-latex*1

2022年4月20日2022年8月20日

模板一

模板二

模板三

第一篇小作文get，十分感谢潘导和任伟师兄，帮助我一个本科生干这些事。阿里嘎多（磕头）

还要继续努力 up up~

附一个文章PDF

Genome-Wide Association Study of Root System Architecture in Maize

2022年1月4日2022年8月20日

【毕业开题】开题研究进度信息详情

该页面作为开题答辩的补充页，对答辩研究进度部分信息进行补充，只展示部分关键数据和代码，更细节部分请问我

下划线文字附带超链接，请自行浏览。只展示关键分析内容和结果，最新进展请以本人的嘴为准。

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-01

1.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_check
```
2.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个候选基因，并对其中高可行度的基因进行连锁分析、单倍体分析、基因表达分析等。

2021年12月28日2022年8月20日

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

2021年12月11日