搜索
楼主: Jimmy

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

  [复制链接]

4

主题

8

帖子

239

积分

中级会员

Rank: 3Rank: 3

积分
239
发表于 2017-1-15 13:02:26 | 显示全部楼层
081-perl+python
最近这周忙着期末考试各种复习没什么时间学习,争取下个星期能把作业补上
回复 支持 反对

使用道具 举报

2

主题

34

帖子

777

积分

高级会员

Rank: 4

积分
777
发表于 2017-1-15 13:04:02 | 显示全部楼层
btrainee 发表于 2017-1-15 10:37
很厉害!大神。
运算速度很快
只是有一点输出的结果里第一条染色体没有统一成大写字符 ...

谢谢提醒,忘记转化第一条染色体的名字了,其他的结果还是正确的
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-1-15 14:52:19 | 显示全部楼层
本帖最后由 anlan 于 2017-1-15 15:01 编辑

034-perl+shell

test.fasta为测试数据

perl代码如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

open my $fh, "test.fasta" or die "Cannot open it";
my ($chr,$seq);
my %hash;
while(<$fh>){
        chomp;
        if ($_ =~ /^>(.*)/){
                $chr = $1;
                $seq = "";
        }else{
                $seq .= $_;
                $hash{$chr} = $seq;
        }
}
close $fh;
map{
        my $len = length($hash{$_});
        my $count_A = $hash{$_} =~ tr/Aa/Aa/;
        my $count_T = $hash{$_} =~ tr/Tt/Tt/;
        my $count_C = $hash{$_} =~ tr/Cc/Cc/;
        my $count_G = $hash{$_} =~ tr/Gg/Gg/;
        my $count_N = $hash{$_} =~ tr/Nn/Nn/;
        my $percent_GC = ($count_C+$count_G)/($count_A+$count_T+$count_C+$count_G);
        print "$_\t$count_A\t$count_T\t$count_C\t$count_G\t$count_N\t$len\t$percent_GC\n";
}sort keys %hash;

perl结果如下:
chr_1        13        10        7        10        4        44        0.425
chr_2        11        6        5        8        4        34        0.433333333333333
chr_3        10        6        3        10        4        33        0.448275862068966
chr_4        9        4        2        7        3        25        0.409090909090909

shell代码如下:

$ awk 'BEGIN{RS=">";FS="\n"} NR>1{seq="";for(i=2;i<=NF;i++) seq=seq""$i;len = length(seq);count_A=0;for(i=0;i<=len;i++){tmp = substr(seq,i+1,1);if(tmp ~ /[Aa]/) count_A++};print $1": "len " "count_A}' test.fasta

结果如下:

chr_1: 44 13
chr_2: 34 11
chr_3: 33 10
chr_4: 25 9

第一列是染色体,第二列碱基长度,第三列是碱基A的数目。
回复 支持 1 反对 0

使用道具 举报

13

主题

33

帖子

243

积分

版主

Rank: 7Rank: 7Rank: 7

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

下载了R和python的视频,想着对R熟悉些,就把R的视频看了下,发现代码等看不清;

于是又来论坛看大家的回复,发现两个有用的提示,一是你的回复,也就尝试了下,发现能成功。
代码如下:
[AppleScript] 纯文本查看 复制代码
setwd("C:/Users/Administrator/Desktop/learn/lecture2/")
rm(list=ls())

source("https://bioconductor.org/biocLite.R")
biocLite("BSgenome.Hsapiens.UCSC.hg19")
biocLite("Biostrings")

library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg19)
alphabetFrequency(Hsapiens$chrY)wd
GC_content<-letterFrequency(Hsapiens$chrY,letters = "CG")
G_content<-letterFrequency(Hsapiens$chrY,letters = "G")
C_content<-letterFrequency(Hsapiens$chrY,letters = "C")
GC_content

GC_pencentage <- letterFrequency(Hsapiens$chrY,letters = "CG")/letterFrequency(Hsapiens$chrY,letters = "ACGT")
GC_pencentage



不过不是很懂 “  alphabetFrequency(Hsapiens$chrY)wd  ” 中wd的具体意思,结合输出数据,整体意思是统计A,G,C,T,N的个数(当然还有些其他字母,不过都是0)。。。可能和wd相关。。。此外试了下chr1,chrx等都没有成功。。。发现删除wd居然可以了,见下图,

123.png

当然除了CG含量,还可以统计其他的含量如A T G C的其他含量

此外,还发现 length(Hsapiens$chr1) 可能统计染色体的长度,这么说R也是很强大的。。。归功于这个包

> length(Hsapiens$chr1)
[1] 249250621


另外一个是,说python的视频很好,也很清楚,于是认真看了下,原理是看懂了,这几天试下。。。

结合上一课的教程,慢慢对Python又些眉目了,体会到了Python在数据处理方面的强大实力了,还是要系统的学习下。

继续跟着生信技能树学习,年底work station 应该会到位,到时候好好整整!!!


回复 支持 反对

使用道具 举报

1

主题

10

帖子

257

积分

中级会员

Rank: 3Rank: 3

积分
257
发表于 2017-1-15 16:25:01 | 显示全部楼层
008+python+R-西花厅

学习了,python教程和直播视频,试着写了代码如下

#!/usr/bin/env python

import re
import sys
import copy
import os

from collections import OrderedDict

#os.chdir("C:\\Users\\a105\\Desktop")

input_file_name="sample.txt"

file=open(input_file_name,"r")

temstring=""
chr_dict=OrderedDict()

for line in file:
        line=line.strip()
        if line.startswith(">"):
                temstring=line
                chr_dict[temstring]=""
        else:
                chr_dict[temstring]+=line
#print(chr_dict)

for seqName,seq in chr_dict.items():
        seqlen=len(seq)
        N_ratio=seq.count("N")/seqlen
        GC_number=seq.upper().count("G")+seq.upper().count("C")
        GC_ration=GC_number/(seqlen-seq.count("N"))
        print(seqName,N_ratio,GC_ration)

这个代码可以用于示例的文件。。。
>chr_1 0.09090909090909091 0.425
>chr_2 0.11764705882352941 0.43333333333333335
>chr_3 0.12121212121212122 0.4482758620689655
>chr_4 0.12 0.4090909090909091
[Finished in 0.1s]
回复 支持 反对

使用道具 举报

0

主题

6

帖子

337

积分

中级会员

Rank: 3Rank: 3

积分
337
发表于 2017-1-15 16:25:54 | 显示全部楼层
035-python+shell

[Shell] 纯文本查看 复制代码
for i in $(seq 1 2) X M;
do 
   cat chr${i}.fa | awk '$1~/^[A-Za-z]/{print $0, length($0)}' | 
   awk -v n=$i '{sum += $2}{s+=gsub(/[GgCc]/,"&")}{sn+=gsub(/N/,"")};END{print "chr"n, sum, s, sn, s/(sum-sn)}'  
done
回复 支持 反对

使用道具 举报

0

主题

5

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-1-15 16:35:59 | 显示全部楼层
168-python-沙鸥
继续使用jupyter notebook,听了两遍的视频,那个84分钟的视频,讲得搜索方法挺不错的。
计算test的文档,将文档逐行读。
import sys
from collections import OrderedDict
sys.chdir('D:/mybioinfor/biotask')

chr_dict = OrderedDict()  #有序的字典
temp_chr = ''   

with open('h19_test.fasta', 'r') as f:
    for line in f:
        line = line.strip()
        if line.startswith('>'):  #以>开头
            temp_chr = line
            chr_dict[temp_chr] = '' #导入到字典
        else:
            chr_dict[temp_chr] += line

with open('atgc_count.txt', 'at') as outcome:
    print('seqID\tseqLength\tcount_A\tcount_T\tGC_content', file=outcome)
    for seqName, seq in chr_dict.items():
        seqLen = len(seq)
        seqID = seqName[1:]
        A = seq.count('A') + seq.count('a')
        T = seq.count('T') + seq.count('t')
        G = seq.count('G') + seq.count('g')
        C = seq.count('C') + seq.count('c')
        N = seq.count('N')
        GC = seq.count('G') + seq.count('g') + seq.count('C') + seq.count('c')
        print('{0}\t{1}\t{2}\t{3}\t{4}'.format(seqID, seqLen, A, T, '%.3f'%(GC/(seqLen-N))), file=outcome)
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-1-15 17:31:15 | 显示全部楼层
本帖最后由 旭日早升 于 2017-3-11 21:44 编辑

201-python-无名

自己的思考解答:学习了第一题之后,感觉python编程处理生物问题挺有意思的。第二题直接对长达3G的genome操作,有点吓人啊。幸好群主给了一个思路,先用小一点的文件测试下,输出结果为A、T、G、C、N的个数。我的思路:初始化一个字典,逐行读入文件并判断,如果是注释行,将染色体名加入字典并初始化一个空字符串,将注释行下的序列赋值给该染色体,文件读完后,对字典进行统计。用测试的小文件还是ok的,但是对基因组操作就一直卡在哪里。先把我的代码放在这里,学习老师的代码后再好好总结下。
我的代码如下:
[Python] 纯文本查看 复制代码
import re
import os 

os.chdir("./")

CHR = {}

with open("UCSC_hg19_chrAll.fasta","rt") as f:
        for line in f:
                if line.startswith(">"):
                        chr_name = line.rstrip().split(">")[1]
                        if chr_name not in CHR:
                                CHR[chr_name] = ""
                else:
                        CHR[chr_name] += line.rstrip()

for key in CHR:
        print(key,len(re.findall("A|a",CHR[key])),len(re.findall("T|t",CHR[key])),len(re.findall("C|c",CHR[key])),len(re.findall("G|g",CHR[key])),len(re.findall("N|n",CHR[key])))
测试文件结果
chr_2 11 6 5 8 4
chr_1 13 10 7 10 4
chr_3 10 6 3 10 4
chr_4 9 4 2 7 3

hg19的结果没有跑出来


学习Teacher Li的视频:李老师讲解了两种方法。其中之一与我的方法类似,但是点拨了一下我们,即不能把所有的序列存入字典,这样占内存很大,可以逐个序列打印,我理解的逐个序列是按染色体来打印,因而结合老师的代码调整了一下,调整代码如下:

[Python] 纯文本查看 复制代码
from collections import OrderedDict

chr_dict = OrderedDict()
temp_chr = ""

with open("UCSC_hg19_chrAll.fasta","rt") as f:
        for line in f:
                line = line.strip()
                if line.startswith(">"):
                        temp_chr = line
                        chr_dict[temp_chr] = ""
                else:
                        chr_dict[temp_chr] += line

for seqName,seq in chr_dict.items():
        A = seq.count("A") + seq.count("a")
        T = seq.count("T") + seq.count("t")
        C = seq.count("C") + seq.count("c")
        G = seq.count("G") + seq.count("g")
        N = seq.count("N") + seq.count("n")
        print(seqName,A,T,C,G,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

但是hg19的文件同样还没有运行出来。因而老师提供了第2种方法,使用pysam包。
由于服务器默认了python2.7,但我一般使用python3,但是import pysam的时候总是报错“ImportError: dynamic module does not define module export function (PyInit_libchtslib)”,一直未找到解决方法,求大神指教。本来想在python2.7下运行的代码,发现pysam的版本太老,没有老师讲的那些references和lengths的命令,我还没有权限,惆怅啊。继续寻找方法解决。
此外,老师在计算CG含量的时候提醒去除N的比例,因为N是gap,所以不能把N计算在总长度内,get到一个小知识点。


终于解决了pysam的问题,这历程不可谓不艰辛啊。首先说明一下报错的问题,询问大神,因为是python2.7的包自动加载到环境变量,所以即使打开了python3也会导入python2.7的模块路径,解决方法是删除python路径(export PYTHONPATH=“”)。这时打开python3后再导入pysam会提示没有该模块。因而尝试通过pip3安装pysam,这时候有报错了,gcc的问题 ,继续询问大神,htslib的C库问题,需要自己导入。接着安装又遇到一个坑,samtools error,在网上查了下需要在samtools的安装路径下建立一个配置文件,我去我没有权限。解决方法,我自己安装了samtools并且优先调用我自己的samtools。接着安装,还是权限问题,这下没办法了,联系python拥有着,搞定。虽然三言两语写出了这历程,但是其间真的不知尝试多少次,还尝试了用source code的方式安装,最终结果还是喜人的。下面就来实践老师的方法吧。
使用pysam的代码
[Python] 纯文本查看 复制代码
import pysam
import os 

os.chdir("./")

hg19 = pysam.FastaFile("UCSC_hg19_chrAll.fasta") 
dir(hg19) #list methods
list(zip(hg19.references,hg19.lengths))
for seqName in hg19.references:
        seq = hg19.fetch(seqName)
        A = seq.count("A") + seq.count("a")
        T = seq.count("T") + seq.count("t")
        C = seq.count("C") + seq.count("c")
        G = seq.count("G") + seq.count("g")
        N = seq.count("N") + seq.count("n")
        print(seqName,A,T,C,G,N)        


运行hg19全基因组的结果:
chr1 65570891 65668756 47024412 47016562 23970000
chr2 71102632 71239379 47915465 47947042 4994855
chr3 58713343 58760485 38653197 38670110 3225295

......
使用pysam真是快啊,简直不是一个等级的。



#20170119#又学习了TeacherDong的86分钟视频,对于编程新手的我来说收获颇多。授人以鱼不如授人以渔,老师这一讲正是给很多从1开始写python的人指明了方法。我学到了什么:首先是一个好习惯,老师编写程序之前写了一个解析问题的大纲,把问题说明、文件格式、数据结构、效率方面等都一一写明,这不仅对于后来编程很有帮助,好习惯对于追踪问题解析问题都很有帮助。其次是利用资源。利用论坛资源搜索功能(例如如何下载hg19),利用网络资源检索命令(例如如何打开文件)。如果遇到一个基本问题或者需要查询一个编程命令,你可以搜索论坛或者网络查找学习或借鉴。另外需要注意生物人编程问题。由于不是计算机科班出身,对于很多问题可能有所忽略。老师提到了代码重复问题以及效率问题(高阶)。代码重复问题是我的一个普遍问题,例如计算A、T、C、G、N的含量,每次都要重复代码,如果使用数组把ATCGN存起来,每次遍历一次数组,使用相同的代码即可。另外效率问题是很多初学者没法照顾的,例如读文件,可以通过readline读入文件,但是怎样快速读取并统计大的文件则需要有一定积累才能解决,老师使用了buffer的方法(搜索快速统计行数找到)。最后是进阶之路。老师给这个课程命了个名字“从1开始写python”。0.5 开始课程,1 基础语法,2 编程技巧,3 效率。开始课程前需要一定的基础,第一层次要有基本语法,第二层次讲究编程技巧,第三层次要提高效率。从现在开始一步步向高阶进发吧。
最后写上老师留的思考作业,用递归的方式解决代码对于小染色体的bug,运行了一下可以,但是感觉有点绕,先把递归结果输出为tuple,然后再转list,还没想到怎么把它直接弄成list。
代码如下:
[Python] 纯文本查看 复制代码
import sys
import time 
import re

start = time.clock()

args = sys.argv

sum_atcg = {}

bases = ["A", "T", "C", "G", "N"]

def get_chr(buffer):
        (buffer,tmp) = buffer.split(">",1)
        if ">" not in tmp:
                return (buffer, tmp)
        else:
                get_chr(tmp)
                return (buffer, get_chr(tmp))

def get_list(tmp):
        tmp_list = []
        while len(tmp) > 1 and type(tmp) == tuple:
                tmp_list.append(tmp[0])
                tmp = tmp[1]
        tmp_list.append(tmp)
        return tmp_list
                
with open(args[1], "r") as Fin:
        tmp = Fin.readline()
        chr_id = re.split(r"\s", tmp)[0][1:]
        sum_atcg[chr_id] = {}
        for base in bases:
                sum_atcg[chr_id][base] = 0
        while 1:
                buffer = Fin.read(1024 * 1024)
                if not buffer:
                        break

                if ">" in buffer:
                        #print(get_chr(buffer))
                        (buffer ,tmps) = get_chr(buffer)
                        #print(tmps,len(tmps))
                        if len(tmps) > 1:
                                tmps = get_list(tmps)
                        #print(tmps)
                        buffer = buffer.upper()
                        for base in bases:
                                sum_atcg[chr_id][base] += buffer.count(base)
                        for tmp in tmps:
                                (tmp, buffer) = tmp.split("\n", 1)
                                chr_id = re.split(r"\s", tmp)[0]
                                sum_atcg[chr_id] = {}
                                for base in bases:
                                        sum_atcg[chr_id][base] = 0
                                buffer = buffer.upper()
                                for base in bases:
                                        sum_atcg[chr_id][base] += buffer.count(base)
                else:
                        buffer = buffer.upper()
                        for base in bases:
                                sum_atcg[chr_id][base] += buffer.count(base)

for chr_id, atcg_count in sum_atcg.items():
        GC = atcg_count["G"] + atcg_count["C"]
        SUM = sum(atcg_count.values())
        print(chr_id)
        for base in bases:
                print("%s: %s" % (base, atcg_count[base]))
        print("SUM:  %s" % (SUM))
        print("GC: %s" % (GC/SUM))
        print("N: %s" % (atcg_count["N"]/SUM))

end = time.clock()

print("used %s s" % str(end - start))

运行结果:
chr8
A: 42767293
T: 42715025
C: 28703983
G: 28702621
N: 3475100
SUM:  146364022
GC: 0.3922180001312071
N: 0.023742856697392477
。。。
used 33.95 s

花费时间还挺少的。


课后思考:TeacherLi提醒了我们要善于利用已经存在的包,不仅是避免重复造轮子,而且效率奇高啊,当然新手练习需要重基础开始。另外,TeacherDong同样从效率出发,告诉了我们,每一个步骤都可能有高效率的方法,这对于编写程序是很重要的进阶需要。










回复 支持 4 反对 0

使用道具 举报

0

主题

2

帖子

92

积分

注册会员

Rank: 2

积分
92
发表于 2017-1-15 19:22:22 | 显示全部楼层
看了视频,和前面各位大神的答案,简直开眼界。但感觉自己有点思维受限,单独用python写不出新的思路来,自己再想想,1月26号前补完答案。 046 R + python  昶
回复 支持 反对

使用道具 举报

0

主题

5

帖子

247

积分

中级会员

Rank: 3Rank: 3

积分
247
发表于 2017-1-15 19:50:10 | 显示全部楼层
064-perl    检查作业的时候,会批注吗?
<#!/usr/bin/perl -w
use strict;

open IN, "hg19.fasta" or die "$!";
open OUT, ">N_GC.txt";

my ($tmp,$keys,@gene,$k,%count,%N_count,%GC_count,%n,%gc);
while (<IN>)
{
        chomp;
        my $line = $_;
        if ($tmp = />chr(\w+)/)
        {
                $keys = $1;
        }
        else
        {
                @gene = split//,$line;
                foreach my $k (@gene)
                {
                        $count{$keys}++;
                        $N_count{$keys}++ if ($k =~ tr/N|n//);
                        $GC_count{$keys}++ if ($k =~ tr/G|C|g|c//);

                }
        }
}
    print OUT "chromosome\tLength\tN_count\tN%\tGC_count\tGC%\n";
        foreach my $keys (sort {$a<=>$b} keys %count)
        {
                 $n{$keys} = sprintf ("%.2f",100*$N_count{$keys}/$count{$keys});
                 $gc{$keys} = sprintf ("%.2f",100*$GC_count{$keys}/$count{$keys});
                print OUT "chr$keys\t$count{$keys}\t$N_count{$keys}\t$n{$keys}\t$GC_coun
t{$keys}\t$gc{$keys}\n";
        }

        close IN;
        close OUT;
>
结果:chromosome      Length  N_count N%      GC_count        GC%
chrX    155270560       4170000 2.69    59679184        38.44
chrM    16571           0.00    7372    44.49
chrY    59373566        33720000        56.79   10252459        17.27
chr1    249250621       23970000        9.62    94040974        37.73
chr2    243199373       4994855 2.05    95862507        39.42
chr3    198022430       3225295 1.63    77323307        39.05
chr4    191154276       3492600 1.83    71776628        37.55
chr5    180915260       3220000 1.78    70218569        38.81
chr6    171115067       3720001 2.17    66306710        38.75
chr7    159138663       3785000 2.38    63308649        39.78
chr8    146364022       3475100 2.37    57406604        39.22
chr9    141213431       21070000        14.92   49639471        35.15
chr10   135534747       4220009 3.11    54607071        40.29
chr11   135006516       3877000 2.87    54504836        40.37
chr12   133851895       3370502 2.52    53252045        39.78
chr13   115169878       19580000        17.00   36827474        31.98
chr14   107349540       19060000        17.76   36099079        33.63
chr15   102531392       20836626        20.32   34475969        33.62
chr16   90354753        11470000        12.69   35332028        39.10
chr17   81195210        3400000 4.19    35428296        43.63
chr18   78077248        3420019 4.38    29702356        38.04
chr19   59128983        3320000 5.61    26989400        45.64
chr20   63025520        3520000 5.59    26257240        41.66
chr21   48129895        13023253        27.06   14334933        29.78
chr22   51304566        16410021        31.99   16745219        32.64
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-12-6 07:15 , Processed in 0.131276 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.