【客户福利】获取ASprofile的hdrs文件

ASprofile 是分析可变剪切的软件之一，软件分析需要一个hdrs文件。该软件自带了hg19的hdrs，但是其它版本或者物种的hdrs文件需要自己获取,但软件并未提供该功能脚本。

hg19.fa.hdrs 部分内容如下：

• 染色体名称，序列长度，除去未知碱基序列的长度和物种。

$ head ASprofile.b-1.0.4/hg19.fa.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)
>chr4 /len=191154276 /nonNlen=187661676 /org=Homo_Sapiens(hg19)
>chr5 /len=180915260 /nonNlen=177695260 /org=Homo_Sapiens(hg19)
>chr6 /len=171115067 /nonNlen=167395066 /org=Homo_Sapiens(hg19)
>chr7 /len=159138663 /nonNlen=155353663 /org=Homo_Sapiens(hg19)
>chr8 /len=146364022 /nonNlen=142888922 /org=Homo_Sapiens(hg19)
>chr9 /len=141213431 /nonNlen=120143431 /org=Homo_Sapiens(hg19)
>chr10 /len=135534747 /nonNlen=131314738 /org=Homo_Sapiens(hg19)

网上找到过一个公开的Python脚本，但是它统计的结果有问题，没办法，自己写了。大家如果也找到了那个脚本，注意验证结果可靠性。

本文分享自己写的perl脚本用法。

该脚本不公开，有需要的客户付费达到一定金额后可免费获取。如果仅需要其它基因组版本或着其它物种的hdrs文件，可以付费获取相应的hdrs文件。

特点

• 参数1：输入文件为基因组文件，不支持压缩的fasta文件。

• 参数2：定义结果文件

• 参数3：物种名称，如果含有特殊字符，如括号()，要用单引号‘ 引起来。

• 对于未知碱基的检测，使用的是N与n。

使用

$ perl get_hdrs_zheng.pl GRCh38.primary_assembly.genome.fa Homo_Sapiens_GRCh38.hdrs 'Homo_Sapiens(hg38)'

查看结果文件

$ head -n 10 Homo_Sapiens_GRCh38.hdrs
>chr1 /len=248956422 /nonNlen=230481012 /org=Homo_Sapiens(hg38)
>chr2 /len=242193529 /nonNlen=240548228 /org=Homo_Sapiens(hg38)
>chr3 /len=198295559 /nonNlen=198100135 /org=Homo_Sapiens(hg38)
>chr4 /len=190214555 /nonNlen=189752667 /org=Homo_Sapiens(hg38)
>chr5 /len=181538259 /nonNlen=181265378 /org=Homo_Sapiens(hg38)
>chr6 /len=170805979 /nonNlen=170078522 /org=Homo_Sapiens(hg38)
>chr7 /len=159345973 /nonNlen=158970131 /org=Homo_Sapiens(hg38)
>chr8 /len=145138636 /nonNlen=144768136 /org=Homo_Sapiens(hg38)
>chr9 /len=138394717 /nonNlen=121790550 /org=Homo_Sapiens(hg38)
>chr10 /len=133797422 /nonNlen=133262962 /org=Homo_Sapiens(hg38)

小结

• 很久没试过网上的脚本拿来直接用了，结果不出所料地又是不能用。
• 本来想写R脚本，发现自己写的读入fasta文件函数太吃内存了，笔记本扛不住。。。