注意：最新版的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:
> str(gwasResults)
'data.frame':16470 obs. of5 variables:
$ SNP: chr"rs1" "rs2" "rs3" "rs4" ...
$ CHR: int1 1 1 1 1 1 1 1 1 1 ...
$ BP: int1 2 3 4 5 6 7 8 9 10 ...
$ P: num0.915 0.937 0.286 0.83 0.642 ...
$ zscore: num0.107 0.0789 1.0666 0.2141 0.4653 ...复制代码
How many SNPs on each chromosome?
> as.data.frame(table(gwasResults$CHR))
Var1 Freq
11 1500
22 1191
33 1040
44945
55877
66825
77784
88750
99721
1010696
1111674
1212655
1313638
1414622
1515608
1616595
1717583
1818572
1919562
2020553
2121544
2222535复制代码Creating manhattan plots
Now, let's make a basic manhattan plot.
manhattan(gwasResults)复制代码
文章图片
2016-4-30 16:50:15 上传
下载附件 (145.41 KB)
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:
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"))复制代码
文章图片
2016-4-30 16:51:05 上传
下载附件 (155.07 KB)
Now, let's look at a single chromosome:
manhattan(subset(gwasResults, CHR == 1))复制代码
文章图片
2016-4-30 16:51:06 上传
下载附件 (81.5 KB)
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.
> str(snpsOfInterest)
chr [1:100] "rs3001" "rs3002" "rs3003" "rs3004" "rs3005" ...
> manhattan(gwasResults, highlight = snpsOfInterest)复制代码
文章图片
2016-4-30 16:51:08 上传
下载附件 (147.22 KB)
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):
manhattan(subset(gwasResults, CHR == 3), highlight = snpsOfInterest, xlim = c(200, 500), main = "Chr 3")复制代码
文章图片
2016-4-30 16:51:09 上传
下载附件 (41.56 KB)
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).
# Add test statistics
> gwasResults
head(gwasResults)
SNP CHR BPPzscore
1 rs111 0.9148 0.10698
2 rs212 0.9371 0.07895
3 rs313 0.2861 1.06663
4 rs414 0.8304 0.21413
5 rs515 0.6417 0.46526
6 rs616 0.5191 0.64474复制代码
文章图片
2016-4-30 16:51:11 上传
下载附件 (288.6 KB)
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
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.
qq(gwasResults$P)复制代码
文章图片
2016-4-30 16:51:14 上传
下载附件 (33.3 KB)
We can otionally supply many other graphical parameters.
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)复制代码
文章图片
2016-4-30 16:51:16 上传
【postgwas|postgwas r语言_如何用R绘制GWAS研究的Manhattan图及QQ图】下载附件 (39.45 KB)
推荐阅读
- ggplot|R语言ggplot在一张图里同时画散点图和折线图
- R语言入门课|R语言导入数据文件（数据导入、加载、读取）、使用Hmisc包的spss.get函数导入SPSS中的por格式文件
- R语言入门课|R语言使用dbinom函数生成二项分布密度数据、使用plot函数可视化二项分布密度数据（Binomial Distribution）
- R语言入门课|R语言使用＜-操作符创建新的变量、使用两个数据列通过加和创建新的数据列（sum variables to make new featurs in dataframe）
- 量子计算|CDR（线性模型处理量子噪声处理方法）
- 机器学习实战|四、案例（北京二手房价影响因素分析）
- R语言|R语言数学统计函数和描述性统计函数
- 代码文件备份 | 3-10（对arraymatrix数据进行初步筛选与分组差异表达分析）
- 单细胞测序|经验总结 | R语言批量读取目录下的文件然后按照行名对其进行整合成为data.frame