论文
High-quality de novo assembly of the apple genome and methylome dynamics of early fruit development
https://www.nature.com/articles/ng.3886/
苹果基因组ng.3886.pdf
今天的推文我们来试着复现一下论文中的Figure3c
论文中对应的图注
Distribution of sequence identity values between genomic copies and consensus repeats in the GDDH13 assembly (based on 2,198,722 data points). The relative frequencies per percentage of identity of the Helitron, TIR, LTR, LINE, SINE and unclassified TEs (NoCat) are represented in different colors.
论文中的部分数据存储在这个链接
https://iris.angers.inra.fr/gddh13/the-apple-genome-downloads.html
这个TE注释里有identity这个值（这里我不太确定是不是用到的这个值来画图）
文件的格式
写个脚本把SINE LINE 和Helitron的值提取出来
import sys
import re
input_txt = sys.argv[1]
output_txt = sys.argv[2]
pattern01 = sys.argv[3]
regexp01 = re.compile(pattern01)
pattern02 = "ID=\S+;"
regexp02 = re.compile(pattern02)
fr = open(input_txt,'r')
ID_list = []
for line in fr:
if 'TargetDescription' in line and len(regexp01.findall(line)) >= 1:
ID_list.append(regexp02.findall(line)[0].replace("ID=","").replace(";",""))
print(len(ID_list))
fr.close()
fr = open(input_txt,'r')
pattern03 = "Parent=\S+;"
regexp03 = re.compile(pattern03)
pattern04 = "Identity=\S+"
regexp04 = re.compile(pattern04)
fw = open(output_txt,'w')
for line in fr:
if "Parent" in line and regexp03.findall(line)[0].replace("Parent=","").replace(";","") in ID_list:
if len(regexp04.findall(line)) >= 1:
fw.write("%s=%s\n"%(regexp04.findall(line)[0],pattern01))
fw.close()
运行脚本
python appleNG.py GDDH13_1-1_TE.gff3 line_identity.txt LINE
python appleNG.py GDDH13_1-1_TE.gff3 sine_identity.txt SINE
python appleNG.py GDDH13_1-1_TE.gff3 helitron_identity.txt Helitron
输出文件格式
接下来是画图代码
library(tidyverse)
dfsine<-read_delim("sine_identity.txt",
delim = "=",
col_names = FALSE)
dfline<-read_delim("line_identity.txt",
delim = "=",
col_names = FALSE)
dfhelitron<-read_delim("helitron_identity.txt",
delim = "=",
col_names = FALSE)
df<-bind_rows(dfsine,dfline,dfhelitron)
ggplot(data=df,aes(x=X2,stat(density),color=X3))+
geom_freqpoly(binwidth=1,linewidth=3)+
theme_classic()+
scale_x_continuous(expand = expansion(mult = c(0,0)),
limits = c(60,100))+
scale_y_continuous(expand = expansion(mult = c(0,0)),
limits = c(0,0.1))+
labs(y="Frequency",x="Identity")+
scale_color_manual(values = c("#2d2884","#c2a20c","#6497d0"),
name="Element")
和论文中的图并不能完全对应上，不太清楚论文中是怎么来统计这个值的
怎么根据这个identity的值算插入时间暂时还没有搞明白
这个图的峰和binwidth的值设置是有关系，binwidth如果改动，line的第二个峰也会有影响，不太明白这个参数应该怎么设置
论文中对这个图的描述文字
To investigate the evolutionary history of TEs in the apple genome, we plotted the distribution of identity values between genomic copies and their consensus sequences (Fig. 3c). Distributions for all classes of repeats showed a peak at 77% identity. By considering the mutation rate that has been reported for LTR-RTs in plants (1.3 × 10−8 base substitutions per site per year40,41), we estimated the age of those insertions as described by the International Human Genome Sequencing Consortium42. We concluded that the peak at 77% identity corresponded to an insertion age of around 21 million years ago (Mya) (Fig. 3c).
We also noted a second peak, particularly for LINE elements, at 98% identity that corresponded to a TE burst at ~1.6 Mya
TEs also have an important role in structuring genomes. The in-depth TE annotation we performed showed a major TE burst in apple that we estimated to have happened around 21 Mya. This affected all types of TEs, suggesting that the precursor of the modern apple underwent environmental changes with resulting stresses that led to the activation of these TEs50. The observed TE burst corresponds to the Miocene epoch (23 Mya to 5 Mya) and may coincide with two events: the divergence between pear and apple48 and an uplift event occurring at the Tian Shan mountains51, which cover the region where the ancestor of the apple originates from52. We hypothesize that these TE bursts, which presumably must have been very different in the predecessor of pear and apple, have contributed to the diversification, and possibly even speciation, of these plants.
推文记录的是自己的学习笔记，很可能存在错误，请大家批判着看
