哪里可以获得.hdrs文件？

extract-as 简介

$ extract-as -help
extract-as txptgtf genome_hdrs [-r tmap goldref]

该脚本需要两个参数：

1. 转录本对应的 gtf 文件

2. 基因组长度统计文件，后缀为 .hdrs，内容如下：

   >chr1 /len=249250621 /nonNlen=225280621 /org=Homo_Sapiens(hg19)
   >chr2 /len=243199373 /nonNlen=238204518 /org=Homo_Sapiens(hg19)
   >chr3 /len=198022430 /nonNlen=194797135 /org=Homo_Sapiens(hg19)
   ...
   

   每一行代表一条染色体，分别给出：

   • 总长度
   • 去除 N 碱基之后的长度
   • 物种信息

哪里可以获得.hdrs文件？

这里有个来源于网络的 python3.0 脚本可以输出 .hdrs 文件，代码如下：

# -*- coding: utf-8 -*-
"""
version: python 3.0
usage: python get_ASprofile_ref_hdrs.py path/species.genome.fa species
"""

import sys
import re
import fileinput
import pandas
import os.path

if len(sys.argv) < 3:
    sys.exit("python error")

FA = sys.argv[1]
# species = os.path.basename(FA).split(".")[0]
species = sys.argv[2]

dic_chr = {}
temp_chr = ''

for line in fileinput.input(FA):
    line = line.strip()
    
    # re.compile函数根据包含的正则表达式的字符串创建模式对象。可以实现更有效率的匹配。
    pat = re.compile(r'^>')
    
    # match()从首字母开始开始匹配，string如果包含pattern子串，则匹配成功，返回Match对象；
    # 失败则返回None，若要完全匹配，pattern要以$结尾。
    match = pat.match(line)
    if match:
        a = line.split(" ")[0]
        temp_chr = a
        
        dic_chr.setdefault(temp_chr, []).append(0)
        dic_chr.setdefault(temp_chr, []).append(0)
    else:
        dic_chr[temp_chr][0] += len(line)
        dic_chr[temp_chr][1] += (line.count("N") + line.count("n"))
    
fileinput.close()


import pandas as pd
import numpy as np

df = pd.DataFrame.from_dict(dic_chr, orient='index')
df = df.reset_index()

for i in range(0, len(df)):
    df.iloc[i, 1] = "/len=" + str(df.iloc[i, 1])
    df.iloc[i, 2] = "/nonNlen=" + str(df.iloc[i, 2])
    
df["species"] = "/org=" + species

np.savetxt(species + '.fa.hdrs.txt', df.values, fmt='%s')

写在后面

这个脚本是我很久之前搜索发现，然后保存下来的，不知道为什么现在浏览器检索不出来作者了。

看到文章“可变剪切分析（一）详细教程”下方评论里很多人在问，想发出来，希望可以帮助大家。