搜索
查看: 1472|回复: 104

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

  [复制链接]

276

主题

478

帖子

1731

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1731
发表于 2017-1-5 08:15:09 | 显示全部楼层 |阅读模式
hg19每条染色体长度,每条染色体N的含量,GC含量。(fasta文件的探索)
需要学习fasta格式,自己学习,需要了解字符串计数!
具体参见:Hg19基因组的分析还有:http://www.biotrainee.com/thread-563-1-1.html

高级作业:蛋白编码区域的GC含量会比基因组其它区域的高吗?

如果你不了解什么是基因组,什么是基因版本,请先看;基因组各种版本对应关系
用R语言也是可以做的,但是不推荐对整个hg19来做,可以做一个测试fasta文件的这些探究。
[AppleScript] 纯文本查看 复制代码
>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4
ATCGTCaaaGANNAATGANGgggTA

如果你计算机太烂了,对hg19基因组搞不定,就用这个测试数据也行。

对这个测试数据,我的处理如下:
[Perl] 纯文本查看 复制代码
perl -alne '{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//); }}END{print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}" foreach sort keys  %N_count}' test.fa 

如果是完整的perl代码:
[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"  

}

结果如下;
>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基因组进行统计。
hg19基因组下载方法有二十多种,我们统一用下面这个:
[Shell] 纯文本查看 复制代码
nohup wget [url=http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz]http://hgdownload.cse.ucsc.edu/g ... Zips/chromFa.tar.gz[/url] &

tar zvfx chromFa.tar.gz

cat *.fa > hg19.fa

rm chr*.fa





上一篇:生信编程直播第五题:根据GTF画基因的多个转录本结构
下一篇:生信编程直播第三题:hg38每条染色体基因,转录本的分布
回复

使用道具 举报

0

主题

4

帖子

30

积分

新手上路

Rank: 1

积分
30
发表于 2017-1-14 12:13:22 | 显示全部楼层
Teanna2017 发表于 2017-1-12 23:27
Python-60

hg19的fa下载了两天都没有下载成功,所以先用测试数据进行了运行,结果如下:

终于下载了hg19的数据,但是里面不同的染色体分成了不同的文件,所以,在读取文件的时候需要读取文件夹中的文件,要import os,

在学习了dongye老师的教学视频时候,修改了代码,没有像之前的代码一样将碱基信息放入字典中,而是识别ATCG直接计数,这样计算结果为:

chr1 249250621
N: 0.1
GC: 0.38
AT: 0.49

chr10 135534747
N: 0.03
GC: 0.4
AT: 0.53

chr11 135006516
N: 0.03
GC: 0.4
AT: 0.52

chr12 133851895
N: 0.03
GC: 0.4
AT: 0.54

chr13 115169878
N: 0.17
GC: 0.32
AT: 0.47

chr14 107349540
N: 0.18
GC: 0.34
AT: 0.45

chr15 102531392
N: 0.2
GC: 0.34
AT: 0.43

chr16 90354753
N: 0.13
GC: 0.39
AT: 0.46

chr17 81195210
N: 0.04
GC: 0.44
AT: 0.5

chr18 78077248
N: 0.04
GC: 0.38
AT: 0.53

chr19 59128983
N: 0.06
GC: 0.46
AT: 0.47

chr2 243199373
N: 0.02
GC: 0.39
AT: 0.54

chr20 63025520
N: 0.06
GC: 0.42
AT: 0.5

chr21 48129895
N: 0.27
GC: 0.3
AT: 0.4

chr22 51304566
N: 0.32
GC: 0.33
AT: 0.34

chr3 198022430
N: 0.02
GC: 0.39
AT: 0.55

chr4 191154276
N: 0.02
GC: 0.38
AT: 0.56

chr5 180915260
N: 0.02
GC: 0.39
AT: 0.55

chr6 171115067
N: 0.02
GC: 0.39
AT: 0.55

chr7 159138663
N: 0.02
GC: 0.4
AT: 0.54

chr8 146364022
N: 0.02
GC: 0.39
AT: 0.54

chr9 141213431
N: 0.15
GC: 0.35
AT: 0.46

chrM 16571
N: 0.0
GC: 0.44
AT: 0.56

chrX 155270560
N: 0.03
GC: 0.38
AT: 0.53

chrY 59373566
N: 0.57
GC: 0.17
AT: 0.23

计时:
real        6m19.964s
user        6m6.884s
sys        0m3.009s

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

filepath = '/Users/sunshine/practice/practice2-hg19/chromFa/'
fileDict = []
filenames = os.listdir(filepath)

for f in filenames:
    if f.startswith('.'): #remove hidden files
        pass
    else:
        fileDict.append(f)

chr_dict = collections.OrderedDict()
bases = ['A','T','C','G','N','a','t','c','g','n']

for f in fileDict:
    f = open(filepath+f,'r')
    for line in f:
        line = line.strip('\n')
        if line.startswith('>'):
            temp = line[1:]
            chr_dict[temp] = {}
            for base in bases:
                chr_dict[temp][base] = 0
        else:
            for base in bases:
                chr_dict[temp][base] += line.count(base)
                
for chr_name, atcgn in chr_dict.items():
    length = sum(atcgn.values())
    N = round((atcgn['N']+atcgn['n'])/length,2)
    GC = round((atcgn['G']+atcgn['C']+atcgn['g']+atcgn['c'])/length,2)
    AT= round((atcgn['A']+atcgn['T']+atcgn['a']+atcgn['c'])/length,2)
    print(chr_name, length)
    print('N: '+ str(N))
    print('GC: '+ str(GC))
    print('AT: '+ str(AT))
    print('')
    
f.close()
回复 支持 3 反对 0

使用道具 举报

0

主题

20

帖子

124

积分

注册会员

Rank: 2

积分
124
发表于 2017-1-11 16:53:35 | 显示全部楼层
最近看书,正好看到Biostrings这个包,发现可以解决这个问题!就试了一下,挺好用!把使用过程贴出来供大家学习!(ps:R语言处于学习阶段,暂时还不会编程。)
[AppleScript] 纯文本查看 复制代码
source("https://bioconductor.org/biocLite.R")
biocLite("BSgenome.Hsapiens.UCSC.hg19")
library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg19)
alphabetFrequency(Hsapiens$chrY)
GC_content<-letterFrequency(Hsapiens$chrY,letters = "CG")
GC_content
GC_pencentage <- letterFrequency(Hsapiens$chrY,letters = "CG")/letterFrequency(Hsapiens$chrY,letters = "ACGT")
GC_pencentage



本帖子中包含更多资源

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

x
人生若只如初见!
回复 支持 3 反对 0

使用道具 举报

4

主题

11

帖子

133

积分

注册会员

Rank: 2

积分
133
QQ
发表于 2017-1-15 17:31:15 | 显示全部楼层
本帖最后由 旭日早升 于 2017-1-19 16:30 编辑

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(">"):
                        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)
                        temp_chr = line
                        chr_dict = OrderedDict()
                        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同样从效率出发,告诉了我们,每一个步骤都可能有高效率的方法,这对于编写程序是很重要的进阶需要。










回复 支持 1 反对 0

使用道具 举报

0

主题

4

帖子

38

积分

新手上路

Rank: 1

积分
38
发表于 2017-1-12 23:03:45 | 显示全部楼层
本帖最后由 lilith098 于 2017-1-15 19:22 编辑

python+137
状态:还没有看视频,根据基础内容和第一题学到的内容尝试解答。

解题思路:
1 下载hg19的数据,因计算机配置原因,只选择chrY.fa 进行统计;
2 删除序列的空格
3 忽略开头字符“>”的行;
4 将所有列输出到一个字符串中;
5 针对字符串完成计数和输出。

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

os.chdir('D:/Programing/study/biotrainee/test2')

seq = str()
with open('chrY.fa',"r") as chrY:
    for line in chrY:
        line = line.strip()     #去除空格
        if line.startswith('>'):
            chr = line
        else:
            seq += line    #将每一行的字符串相加
#seqLen = len(seq)
N = seq.count('N')
GC = seq.count('G') + seq.count('g') + seq.count('C') + seq.count('c')
print (chr, seqLen, N, GC)

输出报错:
TypeError                                 Traceback (most recent call last)<ipython-input-34-a1cd4a2cd5f8> in <module>()----> 1 seqLen = len(seq)      2 N = seq.count('N')      3 GC = seq.count('G') + seq.count('g') + seq.count('C') + seq.count('c')TypeError: 'int' object is not callable
分步输出:



结果:seqLen = len(seq) 此步骤报错,报错内容:TypeError: 'int' object is not callable

求教各位,这个地方报错的原因是什么?我实在没有看出来。

另外,本次作业只针对一个chrY.fa,没有考虑到题目中test.fa的类型。下一步将尝试处理此问题。

20170115
学习了超长的视频,代码如下:
[AppleScript] 纯文本查看 复制代码
#!/usr/bin/python
#Filename:test_1.py

import os
os.chdir('D:/Programing/study/biotrainee/test2')
count_ATCG = {}
Bases = ['a', 't', 'c', 'g', 'n']

#count the base and store in dict. the key is chr, the value is another dict stored the base number respectively
# with open(sys.argv[1], 'r') as fasta:

with open("test.fa", "r") as fasta:
    for line in fasta:
        if line.startswith('>'):
            chr_id = line[1:]
            count_ATCG[chr_id] = {}
            for base in Bases:
                count_ATCG[chr_id][base] = 0
        else:
            line = line.lower()
            for base in Bases:
                count_ATCG[chr_id][base] += line.count(base)

#print the dict
for chr_id, ATCG_count in count_ATCG.items():
    count_sum = sum(ATCG_count.values())
    count_GC = ATCG_count['g'] + ATCG_count['c']
    print(chr_id)
    for base in Bases:
        print("%s: %s" % (base, ATCG_count[base]))

    print("Sum: %s" % count_sum)
    print("N %%: %.2f%%" % (ATCG_count['n']*100.00/count_sum))
    print("GC %%: %.2f%%" % (count_GC*100.00/count_sum))


pycharm运行结果:

chr_1

a: 13
t: 10
c: 7
g: 10
n: 4
Sum: 44
N %: 9.09%
GC %: 38.64%
chr_4

a: 9
t: 4
c: 2
g: 7
n: 3
Sum: 25
N %: 12.00%
GC %: 36.00%
chr_2

a: 11
t: 6
c: 5
g: 8
n: 4
Sum: 34
N %: 11.76%
GC %: 38.24%
chr_3

a: 10
t: 6
c: 3
g: 10
n: 4
Sum: 33
N %: 12.12%
GC %: 39.39%

小结:
1 sys.argv[]输入方式
2 循环语句巩固
3 字典的应用
4 2位浮点百分数的输出

问题:
为什么我的输出不是按照正常的顺序输出?




本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

8

主题

27

帖子

174

积分

注册会员

Rank: 2

积分
174
发表于 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

使用道具 举报

0

主题

5

帖子

63

积分

注册会员

Rank: 2

积分
63
发表于 2017-1-15 02:41:14 | 显示全部楼层
028-Python-Ryu
解题思路:
1.读取fasta文件,fasta文件每个不同的序列开头都有一个“>”号,因此可作为不同染色体的标记
2.统计序列长度使用len(),统计N、G、C个数使用.count()
3.一开始时考虑像李恒的readfq那样,先将序列存入一个变量,最后再统计,但是考虑需要内存以及速度较慢,最后直接每行处理,比较节省内存

完成代码后,再看老师的视频,最后再修改了一下自己的代码,主要是格式化输出以及计算GC含量时取出N的数目再计算。老师介绍的pysam很不错,又学习到新的知识!赞!

明天再做一下统计蛋白编码区的题目。

代码如下:
[Python] 纯文本查看 复制代码
hg19_reference='hg19.fa'
# hg19_reference="test.fa"


with open(hg19_reference,'r')as genome,open("summary.txt",'w')as result:

    result.write("\t".join(("chromosome","length","N number","N percent","GC number","GC percent","\n")))
    summary=["",0.0,0.0,0.0] # chrom, length, #N, #GC
    for line in genome:
        # print(line)
        if line.startswith(">"):
            if summary[0] != "":
                len_seq="%d" % summary[1]
                len_N="%d" % summary[2]
                percent_N="%.02f" % (summary[2]/summary[1])
                len_GC="%d" % summary[3]
                percent_GC="%.02f" % (summary[3]/(summary[1]-summary[2]))
                result.write("\t".join((summary[0],len_seq,len_N,percent_N,len_GC,percent_GC,"\n")))
            print(line)
            summary=["",0.0,0.0,0.0]
            summary[0]=line.strip()
        else:
            summary[1]+=len(line.strip())
            summary[2]+=line.count("N")
            summary[3]+=line.count("G")+line.count("C")+line.count("g")+line.count("c")
    len_seq="%d" % summary[1]
    len_N="%d" % summary[2]
    percent_N="%.02f" % (summary[2]/summary[1])
    len_GC="%d" % summary[3]
    percent_GC="%.02f" % (summary[3]/(summary[1]-summary[2]))
    result.write("\t".join((summary[0],len_seq,len_N,percent_N,len_GC,percent_GC,"\n")))


运行时间在我的笔记本上大约2分钟左右

最后结果如下:
[Plain Text] 纯文本查看 复制代码
chromosome	length	N number	N percent	GC number	GC percent
>chr10	135534747	4220005	0.03	54607071	0.42
>chr11	135006516	3877000	0.03	54504836	0.42
>chr11_gl000202_random	40103	0	0	21899	0.55
>chr12	133851895	3370501	0.03	53252045	0.41
>chr13	115169878	19580000	0.17	36827474	0.39
>chr14	107349540	19060000	0.18	36099079	0.41
>chr15	102531392	20836623	0.2	34475969	0.42
>chr16	90354753	11470000	0.13	35332028	0.45
>chr17_ctg5_hap1	1680828	100000	0.06	718692	0.45
>chr17	81195210	3400000	0.04	35428296	0.46
>chr17_gl000203_random	37498	0	0	12482	0.33
>chr17_gl000204_random	81310	0	0	44286	0.54
>chr17_gl000205_random	174588	0	0	72862	0.42
>chr17_gl000206_random	41001	0	0	21768	0.53
>chr18	78077248	3420015	0.04	29702356	0.4
>chr18_gl000207_random	4262	0	0	1873	0.44
>chr19	59128983	3320000	0.06	26989400	0.48
>chr19_gl000208_random	92689	0	0	34908	0.38
>chr19_gl000209_random	159169	0	0	74008	0.46
>chr1	249250621	23970000	0.1	94040974	0.42
>chr1_gl000191_random	106433	0	0	47198	0.44
>chr1_gl000192_random	547496	0	0	227171	0.41
>chr20	63025520	3520000	0.06	26257240	0.44
>chr21	48129895	13023203	0.27	14334933	0.41
>chr21_gl000210_random	27682	100	0	14977	0.54
>chr22	51304566	16410004	0.32	16745219	0.48
>chr2	243199373	4994851	0.02	95862507	0.4
>chr3	198022430	3225294	0.02	77323307	0.4
>chr4_ctg9_hap1	590426	0	0	214768	0.36
>chr4	191154276	3492600	0.02	71776628	0.38
>chr4_gl000193_random	189789	0	0	81224	0.43
>chr4_gl000194_random	191469	0	0	82827	0.43
>chr5	180915260	3220000	0.02	70218569	0.4
>chr6_apd_hap1	4622290	2301543	0.5	1019557	0.44
>chr6_cox_hap2	4795371	0	0	2142565	0.45
>chr6_dbb_hap3	4610396	406094	0.09	1888816	0.45
>chr6	171115067	3720000	0.02	66306710	0.4
>chr6_mann_hap4	4683263	582522	0.12	1810464	0.44
>chr6_mcf_hap5	4833398	1038487	0.21	1719716	0.45
>chr6_qbl_hap6	4611984	316659	0.07	1936945	0.45
>chr6_ssto_hap7	4928567	755016	0.15	1845164	0.44
>chr7	159138663	3785000	0.02	63308649	0.41
>chr7_gl000195_random	182896	0	0	74370	0.41
>chr8	146364022	3475100	0.02	57406604	0.4
>chr8_gl000196_random	38914	0	0	15429	0.4
>chr8_gl000197_random	37175	100	0	20023	0.54
>chr9	141213431	21070000	0.15	49639471	0.41
>chr9_gl000198_random	90085	0	0	34102	0.38
>chr9_gl000199_random	169874	0	0	64407	0.38
>chr9_gl000200_random	187035	0	0	74716	0.4
>chr9_gl000201_random	36148	0	0	21487	0.59
>chrM	16571	0	0	7372	0.44
>chrUn_gl000211	166566	0	0	64475	0.39
>chrUn_gl000212	186858	0	0	82936	0.44
>chrUn_gl000213	164239	0	0	67177	0.41
>chrUn_gl000214	137718	0	0	57182	0.42
>chrUn_gl000215	172545	0	0	72473	0.42
>chrUn_gl000216	172294	0	0	72319	0.42
>chrUn_gl000217	172149	0	0	64709	0.38
>chrUn_gl000218	161147	0	0	67070	0.42
>chrUn_gl000219	179198	0	0	71609	0.4
>chrUn_gl000220	161802	0	0	78417	0.48
>chrUn_gl000221	155397	0	0	60038	0.39
>chrUn_gl000222	186861	0	0	82170	0.44
>chrUn_gl000223	180455	0	0	77956	0.43
>chrUn_gl000224	179693	0	0	77785	0.43
>chrUn_gl000225	211173	0	0	100631	0.48
>chrUn_gl000226	15008	0	0	5857	0.39
>chrUn_gl000227	128374	0	0	52646	0.41
>chrUn_gl000228	129120	0	0	69641	0.54
>chrUn_gl000229	19913	0	0	9986	0.5
>chrUn_gl000230	43691	0	0	18224	0.42
>chrUn_gl000231	27386	0	0	12232	0.45
>chrUn_gl000232	40652	0	0	17002	0.42
>chrUn_gl000233	45941	0	0	19489	0.42
>chrUn_gl000234	40531	0	0	17457	0.43
>chrUn_gl000235	34474	0	0	13099	0.38
>chrUn_gl000236	41934	0	0	17458	0.42
>chrUn_gl000237	45867	0	0	21403	0.47
>chrUn_gl000238	39939	0	0	15976	0.4
>chrUn_gl000239	33824	0	0	15357	0.45
>chrUn_gl000240	41933	0	0	17842	0.43
>chrUn_gl000241	42152	0	0	15717	0.37
>chrUn_gl000242	43523	0	0	21063	0.48
>chrUn_gl000243	43341	0	0	19940	0.46
>chrUn_gl000244	39929	0	0	17421	0.44
>chrUn_gl000245	36651	0	0	13309	0.36
>chrUn_gl000246	38154	0	0	14746	0.39
>chrUn_gl000247	36422	0	0	15880	0.44
>chrUn_gl000248	39786	0	0	18165	0.46
>chrUn_gl000249	38502	0	0	18011	0.47
>chrX	155270560	4170000	0.03	59679184	0.39
>chrY	59373566	33720000	0.57	10252459	0.4


代码和结果已上传Github,欢迎小伙伴们与我讨论,谢谢。

本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

0

主题

4

帖子

40

积分

新手上路

Rank: 1

积分
40
发表于 2017-1-14 20:47:10 | 显示全部楼层
003+python
我觉得这道题的难点在于如何存储序列名称和序列,并将它们对应起来。
本人还是python菜鸟,还没有看完廖雪峰老师的网站。只是将视频好好看了几遍,根据第一节课老师的推荐,使用的jupyter notebook编辑的,个人感觉确实很好用,支持Tab键自动补齐,自动语法高亮。与第一节课相比,对for循环,if语句,模块等概念和用法有了较为清晰的认识,对于有序字典这个概念,虽然有了认识,但是对于用法还是不太熟悉。需要进一步学习。

对于老师讲的pysam模块,让我进一步了解到了Python的便捷之处,目前已在虚拟机的linux上安装成功,待进一步学习。

本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

0

主题

2

帖子

49

积分

版主

Rank: 7Rank: 7Rank: 7

积分
49
发表于 2017-1-12 22:52:19 | 显示全部楼层
写了一个脚本,输入fasta文件,输出'count_atgc.txt'用来计算A,T,G,C的含量和序列长度
[Python] 纯文本查看 复制代码
#!/usr/bin/env python3
#
# Usage example:
#     count_atgc.py tmp.fasta
# 
# Then, you can open the file "atgc_count.txt"

import sys
args = sys.argv

from collections import OrderedDict

chr_dict = OrderedDict()
temp_chr = '' 

with open(args[1], '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)


得到的文本结果如下图


本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

1

主题

38

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-11 17:19:58 | 显示全部楼层
x2yline 发表于 2017-1-10 23:59
077 R-python-shell-x2yline
(1)R语言实现

要买个好电脑,我的电脑写的程序很多都跑不动
回复 支持 1 反对 0

使用道具 举报

1

主题

12

帖子

114

积分

注册会员

Rank: 2

积分
114
发表于 2017-1-5 21:25:21 | 显示全部楼层
perl:043
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl

open IN,"<$ARGV[0]";
$/ = ">";<IN>;
print "chr\tN(%)\tGC(%)\n";
while(<IN>){
        s/>//;
        my($chr,$seq) = split;
        $len = length $seq;
        $n = ($seq =~ s/N/N/g);
        $g = ($seq =~ s/G/G/g);
        $c = ($seq =~ s/C/C/g);
        $N = ($n/$len) * 100;
        $GC = (($g+$c)/$len) * 100;
        printf "%s\t%.2f\t%.2f\n","$chr","$N","$GC";
}


本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

4

主题

27

帖子

141

积分

注册会员

Rank: 2

积分
141
发表于 2017-1-6 19:13:45 | 显示全部楼层
本帖最后由 babytong 于 2017-1-8 00:50 编辑

025 perl+python

[Python] 纯文本查看 复制代码
def count_N_CG(seq):    length=len(seq[1])
    print("%s \n N: %f \n C+G:%f "%(seq[0],float(seq[1].count('N'))/length,float(seq[1].count('C')+seq[1].count('G'))/length))
with open ('D:/codinglearning/question2/seq.txt') as f:
    all_seq=[]
    all_name=[] 
    seq=''
    for line in f:
        if line.startswith('>'): 
            name_one=line.strip()
            all_name.append(name_one)
            if len(seq)!=0:all_seq.append(seq)            
            seq=''
            continue
        elif len(line.strip()) == 0:#这里用于跳过空行
            continue
        else:
            seq+=line.strip().upper()
    all_name.append(name_one)
    all_seq.append(seq)
    data=zip(all_name,all_seq)
    for i in data:count_N_CG(i)

测试文本的结果如下

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

1

主题

38

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-7 15:03:20 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-7 15:07 编辑

下载了1号染色体的数据进行一个测试,我的电脑h19估计跑不动,但没想到的是chr1一个文件也跑了很久。
解题思路:
1,提取fasta中的序列
2,计算GC和N的含量
[Python] 纯文本查看 复制代码
#coding=utf-8
#第一步提取染色体1的序列
seq=''
with open('E:/chr1.fa/chr1.fa','r',)as f:
    for line in f:
        line=line.rstrip()
        if line.startswith('>'):
            name=line#将第一行染色体的名称存下来
            continue
        seq+=line#将每一行的的核苷酸字符串连接起来(chr1有4985015行)
#第二部计算GC,N含量
#文件中包含很多小写的a,c,t, g因此全部改为大写形式方便统计
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'))
回复 支持 反对

使用道具 举报

4

主题

27

帖子

141

积分

注册会员

Rank: 2

积分
141
发表于 2017-1-8 00:51:15 | 显示全部楼层
Jimmy 发表于 2017-1-7 11:57
重新来一个,我修改了测试文件!

已更新代码与结果
回复 支持 反对

使用道具 举报

0

主题

7

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-8 22:01:49 | 显示全部楼层
139-python-佘璇
一开始也是想着楼上兄弟的做法,后来想着麻烦,看了下貌似挺简单的:

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

10

主题

32

帖子

139

积分

注册会员

Rank: 2

积分
139
发表于 2017-1-9 15:10:58 | 显示全部楼层
sxuan 发表于 2017-1-8 22:01
139-python-佘璇
一开始也是想着楼上兄弟的做法,后来想着麻烦,看了下貌似挺简单的: ...

哈哈 简单粗暴版,不过一般机器跑不动这份程序还有,你的upper方法用早了,应该在split('>')之后,count之前做
回复 支持 反对

使用道具 举报

0

主题

20

帖子

124

积分

注册会员

Rank: 2

积分
124
发表于 2017-1-9 21:05:45 | 显示全部楼层
本帖最后由 y461650833y 于 2017-1-9 21:07 编辑

最近一直在看perl小骆驼书,看这个题目的时候,顺便也上网找了一些相关资料,不过对哈希那里,有时候还是有点懵!不过还是模仿着写了出来,不过也是来来回回改了好几次!
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
while(<>){
chomp;
if(/>(.+)/){
$key=$1;
}
else {
$hash{$key}.= $_;
}
 
}
 foreach $key (sort keys  %hash){
 $len=length $hash{$key};     
 $countA=$hash{$key}=~s/A/A/ig;##全局替换,返回的是次数!考虑大小写,开始没考虑,传代码的时候想起来了##
 $countT=$hash{$key}=~s/T/T/ig;
 $countC=$hash{$key}=~s/C/C/ig;
 $countG=$hash{$key}=~s/G/G/ig;
$countN=$hash{$key}=~s/N/N/ig;
$GC=($countC+$countG)/( $len-$countN)*100;
print ">".$key."\n";
print "$len\t$countA\t$countC\t$countG\t$countT\t$countN\t$GC\n";

 
 }



本帖子中包含更多资源

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

x
人生若只如初见!
回复 支持 反对

使用道具 举报

0

主题

2

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-1-9 21:29:48 | 显示全部楼层
一直在恶补基础,第一题没赶上,等学完基础后再补,第二题自己先做题,由于对基因数据还不太了解,所以用了测试数据来进行编码:
def tj():
    chr_1='ATCGTCGaaAATGAANccNNttGTAAGGTCTNAAccAAttGggG'
    chr_2='ATCGAATGATCGANNNGccTAAGGTCTNAAAAGG'
    chr_3='ATCGTCGANNNGTAATggGAAGGTCTNAAAAGG'
    chr_4='ATCGTCaaaGANNAATGANGgggTA'

    chr=(chr_1,chr_2,chr_3,chr_4) #创建一个元组,用于循环
    d={} #创建字典
    ysj='ATCGN' #需要搜索的字符
    for j in range(1,5): #遍历元组
        for i in range(1,6): #遍历搜索的字符
            d.setdefault(ysj[i-1]+str(j),chr[j-1].upper().count(ysj[i-1])) #添加健-值

    return sorted(d.items()) #如果不进行排序,那么字典是无序的,可以用sorted函数对字典进行排序
回复 支持 反对

使用道具 举报

0

主题

7

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-9 21:46:52 | 显示全部楼层
dongye 发表于 2017-1-9 15:10
哈哈 简单粗暴版,不过一般机器跑不动这份程序还有,你的upper方法用早了,应该在split('>')之后,count之 ...

这个upper用后面有啥区别吗   
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2017-1-25 01:10 , Processed in 0.264885 second(s), 37 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.