GWAS实战教程之制作PLINK格式的文件（下）

2022-08-21 17:14:12 浏览数 (1)

在上期推送中，我带领大家制作了表型文件和.map文件，今天我们学习一下如何制作.ped文件，关于.ped文件的信息请参见往期内容GWAS实战之制作PLINK格式的文件（上）。

由于我们使用的数据集中最显著的位点是位于1号染色体的DDR2基因上（PMID：29216386），因此为了处理简单，我只挑选1号染色体上的位点用于分析。代码如下：

代码语言：javascript复制
options(stringsAsFactors=F)
library(data.table)
geno <-fread('GSE148812_genotyping_results.txt', header=T) #读取基因型文件
head(geno)
##[1] 242901   1620 说明有242901个位点，1619个样本
annot <- fread('myWES.map', sep='t',header=F) # 读取之前做的.map文件
head(annot)
代码语言：javascript复制
annot <- annot[which(annot$V1 == 1),] # 选择1号染色体上的位点
geno <-geno[which(geno$ID_REF%in%annot$V2),] # 选择基因型文件中位于1号染色体上的位点
dim(geno)
##[1]24675 1620
markerID <- geno$ID_REF # 提取marker ID
annot <- annot[match(markerID,annot$V2),] # 匹配好marker的顺序
fwrite(annot, 'myWES_chr2.map', sep = 't', col.names=F)
sampleID <- colnames(geno)[-1] # 提取样本ID，去掉ID_REF
mygeno<- geno[,2:ncol(geno)] #提取基因型数据
mygeno<- t(as.matrix(mygeno)) # 对基因型数据转置，使得行代表样本，列代表marker
colnames(mygeno) <- markerID # 给新数据添加列名
mygeno <- as.data.frame(mygeno)
mygeno[1:5,1:5]
代码语言：javascript复制
##这里新建一个列表，用于存储基因型数据，便于后续使用lapply函数
genoList =list()
for ( i in 1:ncol(mygeno) ) {
       genoList[[i]]<- mygeno[,i]
       }
length(genoList)
##[1] 24675
A1 <- lapply(genoList,function(x){unlist(strsplit(x,""))[seq(1,2*dim(mygeno)[1],2)]}) # 循环切割基因型数据
A1 <-as.data.frame(matrix(unlist(A1),byrow=F, ncol=ncol(mygeno))) # 将切割好的第一个等位基因变成数据框（要按列排）
A2 <- lapply(genoList,function(x){unlist(strsplit(x,""))[seq(2,2*dim(mygeno)[1],2)]}) # 同上
A2 <- as.data.frame(matrix(unlist(A2),byrow=F,ncol=ncol(mygeno))) # 同上
A1[1:5,1:5]
A2[1:5,1:5]
代码语言：javascript复制
## 下面一部是循环插入等位基因数据，因为A1的第k（k=1，2，3等）列与A2的第k列属于同一位点（SNP）
pedList <- list()
for ( i in 1:ncol(A1) ) {
       pedList[[2*i-1]]<- A1[,i]
       pedList[[2*i]]<- A2[,i]
       }
myped <-as.data.frame(matrix(unlist(pedList), byrow=F, ncol=2*ncol(mygeno))) # 转化为数据框
rownames(myped) <- sampleID # 添加行名
dim(myped)
##[1] 1619 49350
sample_info <-fread('./sample_info.csv')
sample_info$sex <- -9
sample_info$sex[which(sample_info$`gender:ch1`=="Male")]<- 1
sample_info$sex[which(sample_info$`gender:ch1`=="Female")]<- 2
sample_info$smoking <- -9
sample_info$smoking[which(sample_info$`smoking_status:ch1`=="Non-smoker")]<- 1
sample_info$smoking[which(sample_info$`smoking_status:ch1`=="Smoker")]<- 2
sample_info$FID <- sample_info$id
sample_info$IID <- sample_info$id
sample_info$PID <- 0
sample_info$MID <- 0
myfam <-sample_info[,c('FID','IID','PID','MID','sex' ,'smoking')]
myfam <- myfam[match(rownames(myped), myfam$FID),]
myped <- cbind(myfam, myped)
myped[1:5,1:12]

可以看出，myped文件的前六列就是.fam文件的前六列，后面就是基因型信息。

代码语言：javascript复制
fwrite(myped, 'myWES_chr2.ped', sep=' ',col.names=F)

关于.ped文件的制作就讲到这里，后续我会和大家介绍如何用PLINK做GWAS研究。

map 推送

0 人点赞