一文教你绘制Cell文章中SNP的地理分布图

「写在前面」

通过我们组三代人的不懈努力，终于将这份工作以长文的形式近期发表在国际知名杂志**'CELL'「上，感兴趣的可以去看一下，」'Ca2 sensor-mediated ROS scavenging suppresses rice immunity and is exploited by a fungal effector'**。而我作为试验室的一份子，主要参与了SNP的分型及地理分布方面的工作。在绘图的过程中，得到了我的生信指路人--徐洲更师兄的帮助。

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」】。

然后在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,

我是一次在图书馆做生信分享的时候跟巩翔宇有过一面之缘。后来在健身房发现两个人都爱好健身，于是就变成他带我撸铁，我给他布置生信自学任务，以及帮他克服一些瓶颈。除了，绘制了这个SNP分布图，实际上他也能处理大部分的组学分析，他的生信自学笔记都在 https://www.jianshu.com/u/efb75b8800e7 里面分享，感兴趣的可以看看他的学习历程 （FROM 洲更）

赞 (0)

相关推荐