搜索
楼主: Jimmy

生信编程直播第二题-hg19基因组序列的一些探究

  [复制链接]

0

主题

17

帖子

159

积分

注册会员

Rank: 2

积分
159
发表于 2017-1-11 18:17:08 | 显示全部楼层
Jimmy 发表于 2017-1-11 14:37
如果你想用R语言解决这个问题,请试一试biostring这个包

好的,我试一试,多谢群主老大。
回复 支持 反对

使用道具 举报

0

主题

17

帖子

159

积分

注册会员

Rank: 2

积分
159
发表于 2017-1-11 18:17:57 | 显示全部楼层
y461650833y 发表于 2017-1-11 16:53
最近看书,正好看到Biostrings这个包,发现可以解决这个问题!就试了一下,挺好用!把使用过程贴出来供大家 ...

这个包确实方便呀,赞!
回复 支持 反对

使用道具 举报

0

主题

2

帖子

70

积分

注册会员

Rank: 2

积分
70
发表于 2017-1-11 21:47:24 | 显示全部楼层
052-perl+R
才开始接触生信编程,懂得的还很少,只好根据论坛上的答案来敲了,有点抄作业的感觉
思路:计算ATGC四种碱基加上N得出总长度,再计算GC比例
[Perl] 纯文本查看 复制代码
 #!/usr/bin/perl
 while (<>){
 chomp;

  if(/^>/){
  $chr=$_
 }
 else{
 $A_count{$chr}+=($_=~tr/Aa//);
 $T_count{$chr}+=($_=~tr/Tt//);
 $C_count{$chr}+=($_=~tr/Cc//);
 $G_count{$chr}+=($_=~tr/Gg//);
 $N_count{$chr}+=($_=~tr/Nn//);
  }
 }
 foreach (sort keys %N_count){
 $length =$A_count{$_}+$T_count{$-}+$C_count{$-}+$G_count{$-}+$N_count{$-};
 $gc = ($G_count{$-}+$C_count{$-})/($A_count{$_}+$T_count{$-}+$C_count{$-}+$G_count{$-});
 print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}\t$length\t$gc\n"
 
 }

使用测试数据做的,得到结果,后来用hg19中的chr1做是就出现out of memory 了,果然电脑不行。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

81

积分

注册会员

Rank: 2

积分
81
发表于 2017-1-11 22:05:21 | 显示全部楼层
Python-50
[Python] 纯文本查看 复制代码
from collections import OrderedDict
chr1_dict = OrderedDict()

temp_chr1 = ""

with open('chr1_dome.fa','r') as hg19:
    for line in hg19:
        line = line.strip()
        if line.startswith('>'):
            temp_chr1 = line
            chr1_dict[temp_chr1] = ""
        else:
            chr1_dict[temp_chr1] += line

for seqName , seq in chr1_dict.items():
   seq_length = len(seq)
   seqs = seq.upper()
   N = seqs.count('N')
   GC = seqs.count('G') + seqs.count('C')

print('GC含量是:%.2f'%(GC/seq_length),"\n"'序列长度是:%d'%(seq_length),"\n"'N的含量:%.2f'%(N/seq_length))
回复 支持 反对

使用道具 举报

0

主题

10

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-11 22:56:24 | 显示全部楼层

[Python] 纯文本查看 复制代码
import time
start = time.clock()
seq =''
with open('chr1.fa','r',)as f:
    for line in f:
        line = line.rstrip()
        if line.startswith('>'):
            continue
        seq += line 
    print (len(seq))
seqs=seq.upper()
print (name)
print (seqs.count('N')/len(seqs)),
print ((seqs.count('C')+seqs.count('G'))*1.0/(len(seqs)-seqs.count('N')))

end = time.clock()
print("The function run time is : %.03f seconds" %(end-start))
回复 支持 反对

使用道具 举报

0

主题

2

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2017-1-12 11:10:23 | 显示全部楼层
065-Perl+shell+R作图:
思路:1)首先把每个染色体号作为哈希的键,把序列连在一行作为哈希的值。
2)用正则表达式计算每个碱基的长度,通过循环把每个哈希的键和值输出。
关键在于把染色体号和序列分别存入哈希的键和值,以及用正则匹配去计算各个碱基的长度,有点忘了去网上查了下思路
#!usr/bin/perl -w
use strict;
open IN,"/root/cainiao/data/second.txt";
open OUT,">/root/cainiao/scripts/second_out.txt";
my $a;
my %hash;
my $keys;
my $length;
my $tmp;
while(<IN>){
        chomp;
        if(/>(.*)/){
                $a=$1;}
        else {
                $hash{$a}.=$_;
        }
}
foreach $keys (sort keys %hash){
        $tmp=$hash{$keys};
        $length=length($tmp);
        (my $count_A)=$tmp =~ s/A/A/ig;
        (my $count_T)=$tmp=~ s/T/T/ig;
        (my $count_C)=$tmp=~s/C/C/ig;
        (my $count_G)=$tmp=~s/G/G/ig;
        (my $count_N)=$tmp=~s/N/N/ig;
        print OUT "$keys\t$count_A\t$count_T\t$count_C\t$count_G\t$count_N\n";
}
结果如下:
chr_1   13      10      7       10      4
chr_2   11      6       5       8       4
chr_3   10      6       3       10      4
chr_4   9       4       2       7       3
回复 支持 反对

使用道具 举报

0

主题

25

帖子

162

积分

注册会员

Rank: 2

积分
162
发表于 2017-1-12 11:26:30 | 显示全部楼层
x2yline 发表于 2017-1-11 18:17
这个包确实方便呀,赞!


这个包计算结果 与你的代码是一样的
还是你厉害
人生若只如初见!
回复 支持 反对

使用道具 举报

0

主题

5

帖子

44

积分

新手上路

Rank: 1

积分
44
发表于 2017-1-12 16:36:13 | 显示全部楼层
本帖最后由 scyxr 于 2017-1-12 21:01 编辑

056-Python更新:谢谢000-python-dongye 和jimmy老师的指导,还有同学的帮助。原来是浮点数的问题,我用的是Python 2.7版本,本来想尝试直接将N或者GC含量写成浮点数类型,水平有限,发现行不通,于是就想着用浮点数1.0 * N 和1.0*GC,这样就可以解决问题,也不影响数值本身。PS:这个问题想了有好一会儿,也不断调试了一阵子才想出来的,然后再去QQ群里准备再看看聊天记录的,发现dongye老师之前已经给出了直接的答案,N*1.0
不过还是很高兴,虽然是很小的一个问题,但是自己还是把它解决掉了,在大家的帮助下,我会再接再厉的。
总结:下次再出现相同代码却不能得出相同结果的情况,需要一步一步检查下语法,然后再看看是否是Python版本的问题。

根据老师视频,自己照着来一遍,思路上先读取每一个以chr开头的部分,然后将所有chr的碱基放到一起,最后分别统计N和GC的含量,注意大小写的GgCc都要包含在内。
[Python] 纯文本查看 复制代码
from collections import OrderedDict
chr_dict = OrderedDict()
tem_chr = ''

with open('hg19.fasta') as h:
    for line in h:
        line=line.strip()
        if line.startswith('>'):
            tem_chr = line
            chr_dict[tem_chr]= ''
        else:
            chr_dict[tem_chr] += line
for seqName, seq in chr_dict.items():
    seqLen = len(seq)
    N = seq.count('N')
    G = seq.count('G') + seq.count('g')
    C = seq.count('C') + seq.count('c')
    print(seqName,seqLen,N/seqLen,(G+C)/(seqLen-N))


程序可以运行,每一个染色体总的碱基数目都是对的,但是遇到跟第一题同样的问题,就是,输出结果同样为零,
('>chr_1', 44, 0, 0)
('>chr_2', 34, 0, 0)
('>chr_3', 33, 0, 0)
('>chr_4', 25, 0, 0)

即使把小数点考虑进去也是如此
[Python] 纯文本查看 复制代码
print(seqName,seqLen,'%.4f'%(N/seqLen),'%.4f'%((G+C)/(seqLen-N)))


('>chr_1', 44, '0.0000', '0.0000')
('>chr_2', 34, '0.0000', '0.0000')
('>chr_3', 33, '0.0000', '0.0000')
('>chr_4', 25, '0.0000', '0.0000')
实在搞不懂问题出在哪里了?代码跟老师的也是一样的啊  Answer: 答案在于浮点数与整数相除的问题以及Python版本的问题,我用的是Python2.7
后面将整数变为浮点数以后,问题解决了。
[Python] 纯文本查看 复制代码
print(seqName,seqLen,'%.4f'%(1.0*N/seqLen),'%.4f'%(1.0*(G+C)/(seqLen-N)))

('>chr_1', 44, '0.0909', '0.4250')
('>chr_2', 34, '0.1176', '0.4333')
('>chr_3', 33, '0.1212', '0.4483')
('>chr_4', 25, '0.1200', '0.4091')


回复 支持 反对

使用道具 举报

392

主题

697

帖子

2423

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
2423
 楼主| 发表于 2017-1-12 16:56:11 | 显示全部楼层
scyxr 发表于 2017-1-12 16:36
056-Python
根据老师视频,自己照着来一遍,思路上先读取每一个以chr开头的部分,然后将所有chr的碱基放到 ...

注意缩进。Python的缩进是有意义的。第一个continue后面的语句都应该拿到if循环外
回复 支持 反对

使用道具 举报

0

主题

22

帖子

158

积分

注册会员

Rank: 2

积分
158
发表于 2017-1-12 19:27:50 | 显示全部楼层
scyxr 发表于 2017-1-12 16:36
056-Python
根据老师视频,自己照着来一遍,思路上先读取每一个以chr开头的部分,然后将所有chr的碱基放到 ...

你用的是python2.x还是python3.x版本?这两个版本的除法有一点不同。具体看以下说明和试验。



本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-3-2 02:27 , Processed in 0.029689 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.