在R中進(jìn)行整理，pmap格式在下面帖子中有詳細(xì)介紹。第一列為SNP的ID立倍，第二列為chr，第三列為SNP的位置侣滩，第四列開始為每個性狀的SNP的P值帐萎。
9.3 GWAS：關(guān)聯(lián)分析——EMMAX - 簡書 (jianshu.com)
> pmap<- read.table("pmap.txt", header = T)
> head(pmap)
ID Chr position p
1 Chr1__100201597 Chr1 100201597 0.546670
2 Chr1__100201604 Chr1 100201604 0.332930
3 Chr1__100201625 Chr1 100201625 0.095997
4 Chr1__100201634 Chr1 100201634 0.182580
5 Chr1__100201638 Chr1 100201638 0.087536
6 Chr1__100201672 Chr1 100201672 0.611720
閾值：
9.4 GWAS：顯著性閾值確定——GEC - 簡書 (jianshu.com)
閾值 = 顯著性/有效SNP數(shù)
> threshold <- 1/292493
> threshold
[1] 3.418885e-06
顯著SNP篩選
> sig <- subset(pmap, p <= threshold)
> sig
[1] ID Chr position p
<0 行> (或0-長度的row.names)
即在此性狀中，沒有顯著的SNP胜卤。
圖形繪制
> library("CMplot")
> CMplot(pmap, threshold = threshold, amplify = F, memo = "", file = "tiff", plot.type=c("m","q"))
函數(shù)解釋同9.3 GWAS：關(guān)聯(lián)分析——EMMAX - 簡書 (jianshu.com)
Manhattan plot
QQ plot
黑白曼哈頓圖
> CMplot(pmap, plot.type="m", col=c("grey30","grey60"), LOG10=TRUE, threshold=threshold,threshold.lty=c(1,2), threshold.lwd=c(1,1), threshold.col=c("black"), amplify=TRUE,chr.den.col=NULL, signal.col=c("red"), signal.cex=c(1,1),signal.pch=c(19,19),file="tiff",memo="",file.output=TRUE,verbose=TRUE)
引用轉(zhuǎn)載請注明出處疆导，如有錯誤敬請指出。
