使用R语言包RIdeogram展示blast双序列比对结果

2020-09-03 16:15:24 浏览数 (2)

RIdeogram 是用来展示 染色体组型 （idiogram）的一个R语言包，比如展示snp的密度分布；展示某类基因在染色体上的分布等等。

具体介绍可以参考

https://www.plob.org/article/16335.html

发表的论文是

RIdeogram: drawing SVG graphics to visualize and map genome-wide data on the idiograms PeerJ 2020

附件中还提供了好几份数据集合代码，非常好的学习素材。

我在看这个包的帮助文档的时候也发现了展示2个或者3个基因组共线性分析的图。

需要准备的数据是两个数据框，

第一个是基因组的基本信息，包括

代码语言：javascript复制
> head(karyotype_dual_comparison)
  Chr Start      End   fill species size  color
1  I      1 23037639 969696   Grape   12 252525
2  II     1 18779884 969696   Grape   12 252525
3 III     1 19341862 969696   Grape   12 252525
4  IV     1 23867706 969696   Grape   12 252525
5   V     1 25021643 969696   Grape   12 252525
6  VI     1 21508407 0ab276   Grape   12 252525

第二个是共线性分析的结果

代码语言：javascript复制
> head(synteny_dual_comparison)
  Species_1  Start_1    End_1 Species_2 Start_2   End_2   fill
1         1 12226377 12267836         2 5900307 5827251 cccccc
2        15  5635667  5667377        17 4459512 4393226 cccccc
3         9  7916366  7945659         3 8618518 8486865 cccccc
4         2  8214553  8242202        18 5964233 6027199 cccccc
5        13  2330522  2356593        14 6224069 6138821 cccccc
6        11 10861038 10886821        10 8099058 8011502 cccccc

了解了基本的数据输入格式，那么我们就可以将blast的比对结果转化成这种形式。

首先是使用blast比对 构建数据库

代码语言：javascript复制
makeblastdb -in mt.fasta -dbtype nucl -out mt

比对

代码语言：javascript复制
blastn -query cp.fasta -db mt -outfmt 6 > output.txt

构造RIdeogram的输入数据

df1.txt文件内容

代码语言：javascript复制
Chr,Start,End,fill,species,size,color
I,1,131478,FF9D1E,chloroplast,12,252525
I,1,444567,FF9D1E,mitochondrion,12,252525

output6.txt文件内容

代码语言：javascript复制
NC_044701 NC_044768 96.76 1603 43 7 54477 56074 375799 374201 0.0 2663
NC_044701 NC_044768 83.69 423 37 19 68192 68587 44003 44420 2e-101 370
NC_044701 NC_044768 85.28 326 38 5 66335 66654 360145 360466 1e-88 327
NC_044701 NC_044768 74.07 891 172 45 102938 103801 332888 333746 1e-83 311
NC_044701 NC_044768 90.54 148 11 3 36180 36324 64144 64291 4e-48 193
NC_044701 NC_044768 95.56 90 3 1 49 138 406657 406745 4e-33 143
NC_044701 NC_044768 96.51 86 2 1 110947 111032 420751 420667 2e-32 141
NC_044701 NC_044768 96.43 84 3 0 31738 31821 384266 384349 5e-32 139
NC_044701 NC_044768 94.94 79 4 0 53913 53991 2926 3004 2e-27 124
NC_044701 NC_044768 96.77 62 2 0 105594 105655 110375 110314 2e-21 104
NC_044701 NC_044768 87.34 79 6 3 88210 88284 240028 240106 2e-16 87.9
NC_044701 NC_044768 100.00 31 0 0 10365 10395 63403 63373 2e-07 58.4
NC_044701 NC_044768 96.97 33 0 1 66299 66330 360093 360125 2e-06 54.7

作图代码

代码语言：javascript复制
df1<-read.csv("df1.txt",stringsAsFactors = F)
df1
df2<-read.csv("output6.txt",header=F,sep="t",stringsAsFactors = F)
df3<-df2[,c(3,7,8,9,10)]
df3$fill<-ifelse(df3$V3>90,"0080cc",
                 ifelse(df3$V3<80,"0ab276","e64e60"))

df3$Species_1<-1
df3$Species_2<-1
head(df3)
df4<-df3%>%
  select(Species_1,V7,V8,Species_2,V9,V10,fill)
colnames(df4)<-colnames(synteny_dual_comparison)
ideogram(karyotype =df1 ,
         synteny=df4,output = "1.svg")
rsvg_pdf("1.svg",'3.pdf')

结果

但是这个图有一个缺点是不能体现出两条序列长度的差异，不知道能不能按长度的比列来显示，如何用代码实现暂时还不知道。 想到一个办法是出图后手动加上刻度线

今天先到这里了

svg data drawing graphics txt

0 人点赞