生信技能树-task1-人类基因组外显子区域长度

使用Perl、Python和R统计人类基因组外显子区域长度

最新推荐文章于 2023-08-19 01:07:33 发布

原创

最新推荐文章于 2023-08-19 01:07:33 发布 · 1.1k 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#生物信息学

这篇博客介绍了如何使用Perl、Python和R三种编程语言，从ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt文件中统计人类基因组外显子区域的总长度。通过读取文件、处理数据并利用各自语言的特性去重，最终得出外显子总长度为35185772。三种语言的执行时间分别为0m39.897s、0m43.192s和6.63437s。

目标：统计人类基因组外显子区域长度

题目数据来源为：ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt

perl

open F,"CCDS.current.txt";
while(<F>){#一行一行读取数据
  next if /^#/;#pass掉带#符号开头的行，即首行
  chomp;#去掉末尾的换行符
  @arr=split /\t/;#以制表符切割读取的每一行数据
  next unless /\[(.*?)\]/;#无外显子坐标的行pass
  @start_end=split /,\s+/,$1;#切割文件的第10列
  foreach(@start_end){
    @loci=split /-/;#以“-”切割每个外显子
    foreach(@loci[0]..@loci[1]){$hash{"$arr[0]:$_"}=1;}#以染色体:坐标符号 格式存入perl哈希中，目的去重
  }
}
close F;

$length=keys %hash;#得到哈希键值长度
print "外显子总长度:$length\n";

输出结果

外显子总长度:35185772

用时 0m39.897s

python3

import re

chr_pos = set()#定义set数据,目的去重
with open("CCDS.current.txt", 'r') as file:
    for line in file:
        if re.match('^#',line): 
            continue #跳过首行
        line = line.split("\t")
        if line[9] != "-":#pass掉第10列为"-