计算反转录转座子插入时间二：提取成对LTRs序列

   上一篇讲了反转录转座子插入时间的计算原理，这篇主要是计算过程。首先使用ltr_finder和ltr_harvest对基因组中的反转录转座子进行预测，之后用LTR_retriever整合预测结果，流程主要参考徐洲更先生的博文，输出结果如下：

图源：https://www.jianshu.com/p/f962d5c40fdf

其中*.fa.pass.list.gff3内包含有所有完整的反转录转座子，内容如下：

每个block表示一个反转录转座子，第一行为整个LTR-RT的区域，4、5行分别是5’LTR和3’LTR，根据这两行位置信息，将两段LTRs序列从fasta文件中提取出来。

def readgff3(gff3file):with open(gff3file, "r") as g3f:f = g3f.readlines()out_list = []small_pair = []for line in f:if line.startswith("#"):passif not line.startswith("chr19"):passelse:lin = line.strip().split("\t")chr_name = lin[0]type = lin[2]start = lin[3]end = lin[4]if type == "repeat_region":out_list.append(small_pair)small_pair = []if type == "long_terminal_repeat":ltr_pair = ">" + chr_name + ":" + start + ".." + endsmall_pair.append(ltr_pair)del out_list[0]out_list.append(small_pair)return out_listpair_ltr = readgff3("14_male.fa.pass.list.gff3")
#print(pair_ltr)def readfasta(inputfasta):with open(inputfasta, "r") as file:sequence = {}out_sequence = {}name_line = []out_name = []f = file.readlines()for line in f:line = line.rstrip()if line.startswith(">"):name_line.append(line)sequence[line] = ""currseqname = lineelse:sequence[currseqname] += linefor item in name_line:if item.startswith(">chr19"):start = int(item.strip().split(":")[1].split("..")[0])#if start < 300000:out_name.append(item)for key in sequence.keys():if key in out_name:out_sequence[key] = sequence[key]return out_sequence, out_namefasta_dic, fasta_name = readfasta("14_male.fa.LTRlib.redundant.fa")
#print(fasta_dic)
#print(len(fasta_dic), len(fasta_name))def findfasta(gfflist, fastadic):pair_dict = {}gff_namedic = {}ltr_fasta = []for key in fastadic.keys():#print(key)name_number = key.split("_")[0]gff_namedic[name_number] = keyfor it in gfflist:starts = int(it[0].split(":")[1].split("..")[0])'''if starts > 300000:passelse:'''fiveltr = gff_namedic[it[0]]five_fasta = fastadic[fiveltr]ltr_fasta.append(fiveltr)ltr_fasta.append(five_fasta)threeltr = gff_namedic[it[1]]pair_dict[fiveltr] = threeltrthree_fasta = fastadic[threeltr]ltr_fasta.append(threeltr)ltr_fasta.append(three_fasta)return ltr_fasta, pair_dictpair_ltr, pair_dictionary = findfasta(pair_ltr,fasta_dic)
print(len(pair_dictionary), len(pair_ltr))def write_out(list, dict, file1, file2):with open(file1, "w") as out_f:for i in list:line = i + "\n"#line = "\t".join(str(a) for a in i) + "\n"#print(type(line))out_f.write(line)out_f.close()with open(file2, "w") as out2_f:for key in dict.keys():line1 = key + "\t"out2_f.write(line1)line2 = dict[key] + "\n"out2_f.write(line2)out2_f.close()
#write_out(l, "ltr_redundant_out.txt")
write_out(pair_ltr, pair_dictionary, "14_extract_all.txt", "14_extract_paral_all.txt")

提取结果如下：

在服务器上可以按行把文件拆成小文件，每个文件中包含一对LTRs，也可以在提取的时候就分开写入文件。

计算核酸演化速率的模型有很多种，我们主要使用Kimura 2-Parameter model：
r=K/2T
其中r是核酸演化速率，在近缘物种中较一致；K是每个位点的平均替代数，T是分化时间。想要进一步理解本公式可以参考北京大学的生物演化课程：https://zh.coursera.org/lecture/shengwu-yanhua/8-1-he-suan-yan-hua-de-su-lu-shang-eRCOU
那么当我们已知近缘物种的演化速率，根据K-2模型算出K，即可推算出两条序列的分化时间。K值的计算使用MEGA，由于我们有多条序列，需要成对比对，可以使用MEGA-cc来批量化处理。
下班，具体处理方式下次再更。

计算反转录转座子插入时间二：提取成对LTRs序列相关推荐

计算反转录转座子插入时间三：MEGA批量化处理
得到了成对的LTRs后,可根据两条LTR序列的不同,根据K-2模型计算K--每个核苷酸位点的平均替代数,使用MEGA进行计算.由于文件众多,需要多次处理,使用MEGA-cc进行处理.安装方法参考htt ...
玉米转座子插入型突变体”五折优惠
"玉米转座子插入型突变体"五折优惠一年一度的双十一来啦!中国农业科学院作物科学研究所重大平台中心隆重推出"玉米转座子插入型突变体"5折优惠活动!即日起至12月 ...
如Java8的LocalDate与Date相互转换、ZonedDateTime等常用操作包含多个使用示例、Java8时区ZoneId的使用方法、Java8时间字符串解析成类
下面将依次介绍 Date转Java8时间类操作 ,Java8时间类LocalDate常用操作(如获得当前日期,两个日期相差多少天,下个星期的日期,下个月第一天等) 解析不同时间字符串成对应的Java8 ...
计算鬼成像学习笔记二：二阶关联函数探究
计算鬼成像学习笔记二:二阶关联函数探究 1 一阶关联函数 2 二阶关联函数 3 二阶关联如何重构物体 4 差分鬼成像关联公式 5 归一化鬼成像关联公式 1 一阶关联函数一阶关联函数是光场的电场强度之 ...
GIF动态图怎么提取成图片？分享三个动态图提取成图片的方法
当我们遇到好看的动态图想让他转成图片使用时该怎么办呢,今天分享三个将动态图提取成每一帧图片的方法. 方法一.使用图像编辑器提取操作方法: 1.打开一个支持 GIF 文件格式的图像编辑器,如 Adob ...
java 比较时间时分的大小_java计算时间差及比较时间大小
java计算时间差及比较时间大小 javaz中对日期时间的处理比较多,代码中列出了3中日期时间计算差值的方法. 比如:现在是2004-03-26 13:31:40 过去是:2004-01-02 11: ...
DATESINPERIOD：计算过去某段时间的指标
在商业场景中,有时会需要计算过去某段时间的指标,以此来预测未来的情况.比如用过去三个月的销售均值预测未来一个月的销售额,用其和未来月份的目标值做对比. 那么怎么计算过去三个月的销售额呢? 我们可以用D ...
php把时间格式转换为时间戳,php如何将时间格式转换成时间戳？
php时间格式转换为时间戳的方法:1.使用mktime()将时间转换为时间戳,语法为"mktime(小时.分钟.秒.月.日.年)":2.使用strtime()将字符串表示的日期转换 ...
oracle插入时间报错,Oracle 插入时间时报错：ORA-01861: 文字与格式字符串不匹配的解决办法...
一.写sql的方式插入到Oracle中往oracle中插入时间 '2007-12-28 10:07:24' 如果直接按照字符串方式,或者,直接使用to_date('2007-12-28 10:07 ...
向mysql中插入时间_Java向mysql中插入时间的方法
ava向MySQL插入当前时间的四种方式和java时间日期格式化的几种方法(案例说明);部分资料参考网络资源 java向MySQL插入当前时间的四种方式第一种:将java.util.Date类型的时 ...

计算反转录转座子插入时间二：提取成对LTRs序列

计算反转录转座子插入时间二：提取成对LTRs序列相关推荐

最新文章

热门文章