转自：https://mp.weixin.qq.com/s/F20F_Mn_8etC7dhVXGm46g
1. RICE SNP的使用
由中国农科院作科所牵头的这项3000份栽培稻的测序工作(「Genomic variation in 3,010 diverse accessions of Asian cultivated rice」)是我们做遗传变异的基础，从2018年到现在，引用量已经达到417次，平均一年有140次的引用，令人叹为观止。为了更好地服务人民群众，有心人就以水稻中测得较好的日本晴作为参考基因组，把这些数据充分利用起来，做成了SNP的网站【https://snp-seek.irri.org/】。
下面简单介绍一下如何使用这个网站。
首先选择【「Genotypes」】。
[图片上传中...(image-25d05a-1634205750205-5)]
然后在locus上输入你的基因号，如rod1的基因号是 LOC_Os06 g03810。
输入基因号后，网站会根据你的基因号，填补染色体信息，等它自动填补完，你就可以点「Search」了。很快就会跳出你的结果，最后你把数据download，这里我选择的是CSV格式的。
2.R绘图部分
2.1 数据准备
我们知道，在正式画图之前有一段数据准备工作，要把数据准备成R包能识别的样子。要画SNP的地理分布图，显而易见，主要使用的数据是SNP及其经纬度。
加载R包并读取数据，查看一下有哪些亚分类。
library(openxlsx)
> library(ggplot2)
> library(maps)
> snp <- read.csv('snp3kvars-LOC_OS06G03810-3594761627300423272.csv',sep = ',',header = T)
> colnames(snp)
#查看一下有哪些亚分类
> table(snp$SUBPOPULATION)
> table(snp$SUBPOPULATION)
admix aro aus ind1A ind1B ind2 ind3 indx japx subtrop temp trop
2 103 76 201 209 205 285 475 615 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)
由于数据库没有提供材料产地的经纬度信息，我们只能根据材料的accession，找到材料的收集地。以材料产地国家的首都作为标志【这也就是图中为什么中国的水稻都在北京的原因🤦】，来粗略看一下不同类型SNP的地理分布。
country <- read.xlsx('3K_country.xlsx',sheet = 2)
> snpc <- merge(snpXG,country,by.x = 'ACCESSION',by.y = 'Genetic_Stock_Accession',all.x = T,sort = F)
剩下的就是机械搬运的活，你需要找到这个国家首都的经纬度，然后填到excel中去。这里我们就直接给出了SNP及其经纬度信息，方便大家进行下一步的操作。
关于如何找到的这个自然变异位点，这点就比较有意思，或多或少有运气的成分。我们这个基因非常小，只有182个氨基酸，而这三千份水稻材料中，只有两个变异位点存在于编码区，其中一个还是同义突变，所以我们的落脚点自然而然就放在了另一个非同义突变身上。同时，这一自然变异位点的确认也得到了黄学辉老师方面的帮助和支持,。
2.2 正式绘图
****读取数据并统计位点的频数。这里解释一下地理上的命名规则，东经正数，西经为负数，北纬正数，南纬负数。总体思路是：先把经纬度及SNP作为一组，统计这一组出现的频数；然后再把经度、纬度和变异位点拆分出来，这样就完成了经纬度和变异位点的频数统计。
loc <- read.table("GT.txt",header =T,fill = T)
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
****在毛坯房上添砖加瓦。****这里通过观察，我们把频数分为4个等级：<50;<50<100;<100<150;>150。大一个等级，点就增加0.5倍。
for ( i in 1:nrow(freq)){
pt_color <- ""
pt_size <- 1
if ( freq[i,'type'] == "G"){ #将G指定为蓝色
pt_color <- "#3399FF"
} else{
pt_color <- "orange" #将T指定为橘色
}
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)
}
最终效果如下
然后把图导入到AI中和其他图进行凭借，就得到了Figure 5H
「参考文献」
1.「Genomic variation in 3,010 diverse accessions of Asian cultivated rice」Nature, 557 (2018), pp. 43-49
2.Mansueto, et al. 「Rice SNP-seek database update: new SNPs, indels, and queries」.Nucl. Acids Res.(2017) 45 (D1): D1075-D1081
3.Ca2+ sensor-mediated ROS scavenging suppresses rice immunity and is exploited by a fungal effector,Cell,2021,,ISSN 0092-8674,
