本文参考了橙子牛奶糖的博客，感谢实验室皮卡丘同志的帮助，转载需注明出处
#仅用于记录，仅供参考
准备：Linux系统（软件：plink, twstats, ）；R语言（包：dplyr, ggplot2, qqman）
1. SNP QC
#不包括imputation
plink2 --bfile filename1 --king-cutoff 0.125 --threads20 --out filename2
#去除亲缘
plink --bfile filename2 --keep SampleID.txt --make-bed --out filename3
#提取应保留样本
plink --bfile filename3 --geno 0.02 --hwe 0.000001 --maf 0.01 --make-bed --out filename4
#去除基因缺失率>0.02、哈迪温伯格平衡P值<0.000001、MAF<0.01的SNP
2. PCA
plink --bfile filename4 --indep-pairwise 50 5 0.2 --out filename5
#提取部分数据进行PCA，加快运算速度，实际差距与总体进行PCA不大
plink --bfile filename4 --extract filename5.prune.in --pca 50 --out filename6
#把上一步得到的结果进行PCA
twstats -t twtable -i filename6.eigenval -o filename7
#查看文件，取连续的p值<0.05的PC（一旦有>0.05的，之后<0.05的PC也不取）
#选取PC的方法还有很多，也可以根据碎石图来确定，本文更偏好此方法
3. Covariance
#选取第二步确定的PC，以及根据背景知识认为应该矫正的变量（如最常见的性别、年龄、身高、体重）。此步可以在R里进行，也可在python、excel中进行，本文仅展示利用R进出处理的步骤
COV=read.table('filename',sep=',',head=T,stringsAsFactors=F)
#读取协变量文件，要求含有FID、IID两列，类似性别的变量应为1/2而不是M/F
SPC=read.table('filename6.eigenvec')
#读取PC文件
library(dplyr)
colnames(SPC)[1:4]=c('FID','IID','PC1','PC2')
#根据需要选取的PC数命名列名
mid=merge(x=SPC,y=COV,by=c('FID','IID'))
#以IID为共同列合并两张表格
last=data.frame(mid[,c(1:8,52:55)])
#数字为所需要的列
write.table(last,file='COV.txt',row.names=F,quote=F,sep='')
#保存在当前路径
4. Phenotype
表型频率作图
(a). 连续型变量（箱线图）
R
pdf('boxplot.pdf')
a=read.table('phenotype_filename',sep='',header=T)
a[a=="-9"]<- NA
dim(a)#得到表型行列数
library(ggplot2)
#这里的代码有点丑，见谅，里面都是个人喜好的参数，可以随意调节
for(iin 3:nrow(a))
{b=colnames(a)[i];
d=paste(b,'Frequency',sep='');
e=data.frame(a[,c(1,i)]);
colnames(e)[2]='Phenotype';
print(ggplot(data=e,aes(x=FID,y=Phenotype))+geom_jitter(shape=20,col=gray(0.7))+geom_boxplot(position=position_dodge(1),outlier.shape=NA,alpha=0.2)+ggtitle(d)+theme(plot.title = element_text(hjust = 0.5)))
}
dev.off()
(b). 分类型变量（频数直方图）
#其他的和箱线图一样
print(ggplot(data=ee,aes(x=factor(Phenotype)))+geom_bar(stat="count")+ggtitle(d)+theme(plot.title= element_text(hjust = 0.5)))
5. GWAS
plink2 --bfile filename4 --pheno phenotype_filename --covar COV.txt --logistic hide-covar --ci 0.95 --covar-variance-standardize --threads 10 --out summary_data1
plink2 --bfile filename4 --pheno phenotype_filename --covar COV.txt --linear hide-covar --ci 0.95 --covar-variance-standardize --threads 10 --out summary_data2
#10为线程
6. 绘制曼哈顿图及Q-Q plot
R
a=read.table('summary_data2.glm.linear')
library(qqman)
png('manhattan_plot_name.png',width=1200,height=800)
manhattan(a,chr='V1',bp='V2',snp='V3',p='V14')
dev.off()
#注意logistic与linear的所需行列可能会有所不同，需查看后再进行操作
#qqplot
png('qqplot_name.png',width=1200,height=800)
p=qq(a$V14,main=paste0('Q-Q plot of PhenotypeX GWAS p-values'))
dev.off()
#膨胀系数
p=a$V14
z = qnorm(p/ 2)
lambda =round(median(z^2, na.rm = TRUE) / 0.454, 3)
write.table(lambda,paste0('inflationfactor.txt'),sep='',quote=F,row.names=F)
转载需注明出处
