搜索
查看: 7472|回复: 7

人类hg19这个参考基因组里面各个染色体长度是多少?

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-8-27 21:13:06 | 显示全部楼层 |阅读模式
这个是最基本的问题,可以参考:http://www.bio-info-trainee.com/425.html
首先下载hg19文件,然后写一个脚本统计每条染色体的长度

希望大家都能自己动手练习!



上一篇:试水贴
下一篇:生信编程直播第一题:人类基因组的外显子区域到底有多...
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

2

主题

8

帖子

532

积分

信息监察员

人见人爱的

Rank: 9Rank: 9Rank: 9

积分
532
发表于 2016-9-23 11:16:56 | 显示全部楼层
嘿嘿,我试着用python写了个简单的版本,练练手。

代码如下:
[Python] 纯文本查看 复制代码
import re
import fileinput

hg19 = "./hg19.fasta"
dic_chr = {}
temp_chr = ''

for line in fileinput.input(hg19):
    line = line.strip()
    pat = re.compile(r'^>')
    match = pat.match(line)
    if match:
        temp_chr = line
        dic_chr[temp_chr] = 0
    else:
        dic_chr[temp_chr] += len(line)
fileinput.close()

for item in dic_chr.items():
    print(item)


用时3分钟,不知道是集群慢,还是代码效率低。
但结果是对的。



回复 支持 1 反对 0

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-10-19 18:59:21 | 显示全部楼层
bioxin 发表于 2016-9-23 11:16
嘿嘿,我试着用python写了个简单的版本,练练手。

代码如下:

一般都是单核运算,跟集群木有关系的。但是我们时间差的有点远
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

5

帖子

135

积分

注册会员

Rank: 2

积分
135
发表于 2016-10-20 00:56:45 | 显示全部楼层
      1 start   23-42-52
      2 >chr1   249250621
      3 >chr2   243199373
      4 >chr3   198022430
      5 >chr4   191154276
      6 >chr5   180915260
      7 >chr6   171115067
      8 >chr7   159138663
      9 >chr8   146364022
     10 >chr9   141213431
     11 >chr10  135534747
     12 >chr11  135006516
     13 >chr12  133851895
     14 >chr13  115169878
     15 >chr14  107349540
     16 >chr15  102531392
     17 >chr16  90354753
     18 >chr17  81195210
     19 >chr18  78077248
     20 >chr19  59128983
     21 >chr20  63025520
     22 >chr21  48129895
     23 >chr22  51304566
     24 >chrX   155270560
     25 >chrY   59373566
     26 >chrM   16571
     27 end     23-43-13

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2018-4-25 15:55:16 | 显示全部楼层
用python的交互界面是可以拿到结果的,但是用服务器运行总告诉我 d[seqname] += line.upper().strip()这一项有问题:
Traceback (most recent call last):
  File "len.py", line 11, in <module>
    d[seqname] += line.upper().strip()
KeyError: ''
是否有人遇到过这种情况,很不解!


from collections import OrderedDict
with open(r"./test.fa") as f :
    d = OrderedDict()
    seqname = ""
    for line in f :
        if line.startswith(">"):
            seqname = line.strip()
            d[seqname] = ""
            continue
        else:
            d[seqname] += line.upper().strip()
    print d

    d_length = OrderedDict()
    for key , value in d.items():
        d_length[key] = len(value)
        print key , "length is ",d_length[key]
回复 支持 反对

使用道具 举报

0

主题

2

帖子

107

积分

注册会员

Rank: 2

积分
107
发表于 2018-8-27 18:45:02 | 显示全部楼层
你好 楼主 我看到你的实例中没有统计出chrY的数据,然后我用拟南芥基因组试了下,也无法统计出最后一个染色体的长度,不知道是为什么?
回复 支持 反对

使用道具 举报

3

主题

20

帖子

336

积分

中级会员

Rank: 3Rank: 3

积分
336
发表于 2019-1-11 16:30:28 | 显示全部楼层
bioxin 发表于 2016-9-23 11:16
嘿嘿,我试着用python写了个简单的版本,练练手。

代码如下:

你的代码统计的内容比楼主的更全面
回复 支持 反对

使用道具 举报

0

主题

1

帖子

29

积分

新手上路

Rank: 1

积分
29
QQ
发表于 2019-3-19 17:53:38 | 显示全部楼层
本帖最后由 雁栖湖学dinger 于 2019-3-19 17:55 编辑

用perl写了个简单的脚本,似乎更快

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use File::Basename;
my ($name,$length);
$name = basename($ARGV[0]);
open IN,"$ARGV[0]" or die "Can't read file: $!";
open OUT,"> $name.count" or die "Can't write file: $!";
chomp($name = <IN>);
while (<IN>){
        chomp;
        if (/^\>/){
                print OUT "$name\t$length\n";
                $name = $_;
                $length = 0;
        }else{
                $length += length($_);
        }
}
print OUT "$name\t$length\n";











本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
芷若,oh 我的芷若
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-21 05:35 , Processed in 0.035118 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.