如何用R绘制GWAS研究的Manhattan图及QQ图

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

原文网址：https://cran.r-project.org/web/p ... ignettes/qqman.html
注意：最新版的qqman程序包是基于R 3.2.5版本开发的，所以使用该程序包之前注意更新R软件到最新版本，同时安装最新版的qqman程序包。虽然以前版本的R也能运行相关程序，但是部分功能显示不够完善。

以下内容因为格式问题将源网页内容稍加调整
Intro to the qqman package

The qqman package includes functions for creating manhattan plots and q-q plots from GWAS results. ThegwasResults data.frame included with the package has simulated results for 16,470 SNPs on 22 chromosomes. Take a look at the data:

1. > str(gwasResults)
   
2. 'data.frame':   16470 obs. of  5 variables:
   
3. $ SNP   : chr  "rs1" "rs2" "rs3" "rs4" ...
   
4. $ CHR   : int  1 1 1 1 1 1 1 1 1 1 ...
   
5. $ BP    : int  1 2 3 4 5 6 7 8 9 10 ...
   
6. $ P     : num  0.915 0.937 0.286 0.83 0.642 ...
   
7. $ zscore: num  0.107 0.0789 1.0666 0.2141 0.4653 ...

复制代码

How many SNPs on each chromosome?

1. > as.data.frame(table(gwasResults$CHR))
   
2.    Var1 Freq
   
3. 1     1 1500
   
4. 2     2 1191
   
5. 3     3 1040
   
6. 4     4  945
   
7. 5     5  877
   
8. 6     6  825
   
9. 7     7  784
   
10. 8     8  750
    
11. 9     9  721
    
12. 10   10  696
    
13. 11   11  674
    
14. 12   12  655
    
15. 13   13  638
    
16. 14   14  622
    
17. 15   15  608
    
18. 16   16  595
    
19. 17   17  583
    
20. 18   18  572
    
21. 19   19  562
    
22. 20   20  553
    
23. 21   21  544
    
24. 22   22  535

复制代码

Creating manhattan plots

Now, let's make a basic manhattan plot.

1. manhattan(gwasResults)

复制代码

We can also pass in other graphical parameters. Let's add a title (main=), increase the y-axis limit (ylim=), reduce the point size to 60% (cex=), and reduce the font size of the axis labels to 90% (cex.axis=). While we're at it, let's change the colors (col=), remove the suggestive and genome-wide significance lines, and supply our own labels for the chromosomes:

1. manhattan(gwasResults, main = "Manhattan Plot", ylim = c(0, 10), cex = 0.6, cex.axis = 0.9, col = c("blue4","orange3"), suggestiveline = F, genomewideline = F, chrlabs = c(1:20, "P", "Q"))

复制代码

Now, let's look at a single chromosome:

1. manhattan(subset(gwasResults, CHR == 1))

复制代码

Let's highlight some SNPs of interest on chromosome 3. The 100 SNPs we're highlighting here are in a character vector called snpsOfInterest. You'll get a warning if you try to highlight SNPs that don't exist.

1. > str(snpsOfInterest)
   
2. chr [1:100] "rs3001" "rs3002" "rs3003" "rs3004" "rs3005" ...
   
3. > manhattan(gwasResults, highlight = snpsOfInterest)

复制代码

We can combine highlighting and limiting to a single chromosome, and use the xlim graphical parameter to zoom in on a region of interest (between position 200-500):

1. manhattan(subset(gwasResults, CHR == 3), highlight = snpsOfInterest, xlim = c(200, 500), main = "Chr 3")

复制代码

Finally, the manhattan function can be used to plot any value, not just p-values. Here, we'll simply call the function passing to the p= argument the name of the column we want to plot instead of the default “P” column. In this example, let's create a test statistic (“zscore”), plot that instead of p-values, change the y-axis label, and remove the default log transformation. We'll also remove the genomewide and suggestive lines because these are only meaningful if you're plotting -log10(p-values).

1. # Add test statistics
   
2. > gwasResults <- transform(gwasResults, zscore = qnorm(P/2, lower.tail = FALSE))
   
3. head(gwasResults)
   
4.   SNP CHR BP      P  zscore
   
5. 1 rs1   1  1 0.9148 0.10698
   
6. 2 rs2   1  2 0.9371 0.07895
   
7. 3 rs3   1  3 0.2861 1.06663
   
8. 4 rs4   1  4 0.8304 0.21413
   
9. 5 rs5   1  5 0.6417 0.46526
   
10. 6 rs6   1  6 0.5191 0.64474

复制代码

A few notes on creating manhattan plots:

• Run str(gwasResults). Notice that the gwasResults data.frame has SNP, chromosome, position, and p-value columns named SNP, CHR, BP, and P. If you're creating a manhattan plot and your column names are different, you'll have to pass the column names to the chr=, bp=, p=, and snp= arguments. See help(manhattan) for details.
• The chromosome column must be numeric. If you have “X,” “Y,” or “MT” chromosomes, you'll need to rename these 23, 24, 25, etc. You can modify the source code (e.g., fix(manhattan)) to change the line designating the axis tick labels (labs <- unique(d$CHR)) to set this to whatever you'd like it to be.
• If you'd like to change the color of the highlight or the suggestive/genomewide lines, you'll need to modify the source code. Search for col="blue", col="red", or col="green3" to modify the suggestive line, genomewide line, and highlight colors, respectively.
  

Creating Q-Q plots

Creating Q-Q plots is straightforward - simply supply a vector of p-values to the qq() function.

1. qq(gwasResults$P)

复制代码

We can otionally supply many other graphical parameters.

1. qq(gwasResults$P, main = "Q-Q plot of GWAS p-values", xlim = c(0, 7), ylim = c(0, 12), pch = 18, col = "blue4", cex =1.5, las = 1)

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

关键词：Manhattan Man hat ATT WAS 如何