R绘制部分
1、数据准备
代码来源及原始数据来源cell文章作者，代码及结果部分根据实际情况稍作修改便于演示。代码文件根据你实际情况进行修改，切记！
##加载R包并读取数据，查看一下有哪些亚分类
library(openxlsx)
library(ggplot2)
library(maps)
#下载的文件名记得改下，每次下载数据的名称都不同
snp <- read.csv('snp3kvars-CHR6-23449621-23649621-4238059175866501290',sep = ',',header = T)
colnames(snp)
#查看一下有哪些亚分类
table(snp$SUBPOPULATION)
admix aro aus ind1A ind1B ind2 ind3 indx
2 103 76 201 209 205 285 475 615
japx subtrop temp trop
83 112 288 372
##根据文章的分类标准，我们只挑选籼稻和粳稻的SNP
attach(snp)
snpXG <- snp[SUBPOPULATION == 'ind1A' | SUBPOPULATION == 'ind1B'
|SUBPOPULATION == 'ind2' | SUBPOPULATION == 'ind3'
|SUBPOPULATION == 'indx' |SUBPOPULATION == 'japx'
|SUBPOPULATION == 'temp' |SUBPOPULATION == 'subtrop'
|SUBPOPULATION == 'trop', ]
detach(snp)
country <- read.xlsx('3K_country.xlsx')
snpc <- merge(snpXG,country,by.x = 'ACCESSION',by.y = 'Genetic_Stock_Accession',all.x = T,sort = F)
#材料的经纬度信息利用作者提供的即可，
或者从之前3K发表的文章可以下载，
华中农业大学谢为博老师发表的PNAS上面有具体信息可以下载！
2、R正式绘图
image.png
检查没问题直接往后进行了
##读取数据并统计位点的频数，经纬度需要自己进行转换注意。
loc <- read_excel("GT.xlsx")
info <- paste(loc$lat,loc$lon,loc$SNP, sep = "_")
freq <- as.data.frame(table(info), stringsAsFactors = FALSE)
info_split <- do.call(rbind, strsplit(freq$info, split = "_", fixed = TRUE))
freq$lat <- as.integer(info_split[,1])
freq$lon <- as.integer(info_split[,2])
freq$type <- info_split[,3]
##草图
map('world', fill = TRUE, col = 'grey',resolution=1,
mar = c(0.01, 0.01, par("mar")[2], 0.01),border=NA,
xlim = c(-180,180),ylim = c(-90,90))
map.axes(cex.axis=0.8) #给草图加上axis
image.png
3、图片进行美化
对不同水稻品种ROD1等位基因频率的分析揭示了亚种特有的分布(图S6J)。在分析的803个粳稻品种中，只有5.6%(45/803)携带SNP1A等位基因；相比之下，在1,531个(876/1531)个水稻品种中，该等位基因存在于57.2%(图5C)。为了追踪SNP1的进化起源，我们比较了44份野生稻材料的ROD1序列，发现只有4份材料含有SNP1A，而其他40份材料都含有SNP1C。因此，SNP1A主要存在于水稻中，而SNP1C则主要存在于野生稻和水稻中。SNP1A和SNP1C在两个亚种间分布的显著差异促使我们调查了它们的地理分布。我们发现，含有SNP1A的材料主要集中在热带和亚热带(低纬)地区，而含有SNP1C的材料倾向于在温带地区丰富(图5H)。综上所述，这些结果表明，ROD1的等位基因变异调节了ROS的动态平衡，从而有助于亚种间的抗病。
image.png
##**这里通过观察，我们把频数分为4个等级：<50;<50<100;<100<150;>150。大一个等级，点就增加0.5倍。
#代码部分说明此图只用了原始数据的100份，所以得到的结果有差异。
for ( i in 1:nrow(freq)){
pt_color <- ""
pt_size <- 1
if ( freq[i,'type'] == "A"){ #将A指定为蓝色
pt_color <- "#3399FF"
} else{
pt_color <- "orange" #将C指定为橘色
}
if (freq[i,'Freq'] < 50){ #对频数进行分级
pt_size <- 1
} else if (freq[i,'Freq'] < 50){
pt_size <- 1.5
} else if(freq[i,'Freq'] < 150){
pt_size <- 2
} else{
pt_size <- 2.5
}
points(x=freq[i,'lon'], y =freq[i,'lat'] , #使用内置的点图进行画图，位置为经纬度，大小取决于频数，不同的变异类型填充不同的颜色。
cex = pt_size, col=pt_color, pch = 19)
}
4、最终结果
流程都非常清楚了，代码都是亲测完整流程跑的，其中报错部分都进行调试，进行部分代码修改！
image.png
最终的图片还可以进一步进行加工美化，使得图片更加美观
image.png
数据来源，大家可以自行阅读：
image.png
今天就先给大家介绍到这里，希望大家的科研能有所帮助！祝您科研顺利快乐！
生信漫谈
生信漫谈
生信漫谈，认识生信，学习生信，跨越生信入门路上的障碍，从而利用生信技术解决科研学习路上的绊脚石！
网友评论