全基因组关联分析（GWAS）实操

数据及代码来源：GWAS_in_R下面所用数据均可在此链接下载，作者为githubwhussain2

1.表型数据的读取

设置工作目录

rm(list=ls())

setwd("G:\\GWAS\\R-for-Plant-Breeding-master\\GWAS_in_R")

读取水稻表型数据

> rice.pheno <-read.table("./datafiles/RiceDiversity_44K_Phenotypes_34traits_PLINK.txt",header=TRUE,stringsAsFactors =FALSE,sep ="\t")

#选择适当的列

>rice.pheno <- rice.pheno[,c(2,3,14)]

#将 NSFTV_ 粘贴到每个基因型名称前

> rice.pheno$NSFTVID <- paste0("NSFTV_",rice.pheno$NSFTVID)

删除不存在数据的项

> rice.pheno <- na.omit(rice.pheno)

检查数据的维度

> dim(rice.pheno)

 读取标记数据

从.ped文件中读取基因型信息

> Geno <- read_ped("./datafiles/sativas413.ped")

此时可能会报错

Error in read_ped("./datafiles/sativas413.ped") :

could not find function "read_ped"没有发现函数read_ped

安装R包BGLR可解决问题

>install.packages('BGLR')

载入包

>library(BGLR）

>Geno<-read_ped("./datafiles/sativas413.ped")

>p=Geno$p

>n=Geno$n

>Geno=Geno$x

从.fam文件中载入基因型和编号

>FAM <- read.table("./datafiles/sativas413.fam")

从.map文件中读取数据

>MAP <- read.table("./datafiles/sativas413.map")

在.map文件中重新编码标记数据

> Geno[Geno ==2]<- NA

> Geno[Geno ==3] <- 2

将标记数据转换为矩阵并转置并检查尺寸

> Geno<-matrix(Geno,nrow = p,ncol = n,byrow = TRUE)

> Geno <-t(Geno)

> dim(Geno)

将行名和列名分配给标记文件

> colnames(Geno)<-MAP$V2

添加存储在第二列中的行名并将 NSFTV_ID_ 粘贴到每个行名

> row.names(Geno)<-paste0("NSFTV_",FAM$V2)

2.数据过滤

2.1 标记数据质控

2.1.1 SNP cal rate

检查 snps 的缺失百分比

>miss_snp<- data.frame(SNPs=apply(Geno,2,function(x){ sum(is.na(x))/nrow(Geno)*100}))

绘制缺失数据的直方图

>hist(miss_snp$SNPs,col = "lightblue",xlab = paste("Missing Percenatge"),ylab = "Proportion of markers",main = NULL)

标记大于 25% 的缺失数据

> mrt<-which(miss_snp>25)

从原始文件中删除这些标记并保存

> Geno<-Geno[,-c(mrt)]

> dim(Geno)

>missing_markers<-subset(miss_snp, SNPs>25)

2.1.2 Genotype Call Rate

检查样本中的缺失百分比

> miss_geno<-data.frame(Genotypes=apply(Geno, 1, name <- function(x) {sum(is.na(x))/ncol(Geno)*100 }))

查看条形图，以便更好地可视化样本中的缺失数据

> hist(miss_geno$Genotypes,col = "lightblue",ylab = "Proportion of genotypes",xlab = "Missing Percenatge",main = NULL)

识别缺失数据超过 15% 的样本

> gen<-which(miss_geno>15)

> Geno<-Geno[-c(gen),]

> dim(Geno)

被丢弃的样本列表

> missing_genotypes<-subset(miss_geno,Genotypes>15)

> head(missing_genotypes)

          Genotypes

NSFTV_72   15.97425

NSFTV_194  16.60705

NSFTV_252  18.38992

NSFTV_260  17.04452

2.1.3检查样本间的杂合性

获得杂合百分比

> htero_geno<-data.frame(htr=apply(Geno,1,function(x){(sum((x==1),na.rm = TRUE)/ncol(Geno)*100)}))

使用条形图和直方图进行可视化

> par(mfrow=c(1,2))

> barplot(htero_geno$htr,xlab = "Genotypes",ylab = "Heterozygosity percentage",col = "darkblue")

> hist(htero_geno$htr,xlab = "Heterozygosity percentage",ylab = "Proportion of genotypes",col = "lightblue",main= NULL)

> htrg<-which(htero_geno>4)

> Geno<-Geno[-c(htrg),]

列出被删除的基因型列表

> htrg<-subset(htero_geno,htr>4)

> dim(Geno)

---

写在最后

明天继续看看根据代码重现。从前看山是山，看水是水，现在看山还是山，看水还是水，没有一点点改变，似乎随着时间的流逝，增长的只是我的年龄。