搜索
查看: 2608|回复: 129

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

  [复制链接]

392

主题

697

帖子

2423

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
2423
发表于 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每条染色体基因,转录本的分布
回复

使用道具 举报

4

主题

13

帖子

169

积分

版主

Rank: 7Rank: 7Rank: 7

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










回复 支持 3 反对 0

使用道具 举报

0

主题

6

帖子

44

积分

新手上路

Rank: 1

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

主题

25

帖子

162

积分

注册会员

Rank: 2

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

使用道具 举报

0

主题

6

帖子

55

积分

注册会员

Rank: 2

积分
55
发表于 2017-1-30 01:16:55 | 显示全部楼层
熊海清 发表于 2017-1-30 01:01
074-python

下载了示例demo数据,用视频中老师提到的方法进行了学习

类似方法得到人基因组结果如下,运行非常快
('chrM', 16571, '0.00', '0.44')
('chr1', 249250621, '0.10', '0.42')
('chr2', 243199373, '0.02', '0.40')
('chr3', 198022430, '0.02', '0.40')
('chr4', 191154276, '0.02', '0.38')
('chr5', 180915260, '0.02', '0.40')
('chr6', 171115067, '0.02', '0.40')
('chr7', 159138663, '0.02', '0.41')
('chr8', 146364022, '0.02', '0.40')
('chr9', 141213431, '0.15', '0.41')
('chr10', 135534747, '0.03', '0.42')
('chr11', 135006516, '0.03', '0.42')
('chr12', 133851895, '0.03', '0.41')
('chr13', 115169878, '0.17', '0.39')
('chr14', 107349540, '0.18', '0.41')
('chr15', 102531392, '0.20', '0.42')
('chr16', 90354753, '0.13', '0.45')
('chr17', 81195210, '0.04', '0.46')
('chr18', 78077248, '0.04', '0.40')
('chr19', 59128983, '0.06', '0.48')
('chr20', 63025520, '0.06', '0.44')
('chr21', 48129895, '0.27', '0.41')
('chr22', 51304566, '0.32', '0.48')
('chrX', 155270560, '0.03', '0.39')
('chrY', 59373566, '0.57', '0.40')
回复 支持 1 反对 0

使用道具 举报

0

主题

5

帖子

48

积分

新手上路

Rank: 1

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

主题

32

帖子

214

积分

中级会员

Rank: 3Rank: 3

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

主题

6

帖子

75

积分

注册会员

Rank: 2

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

帖子

56

积分

注册会员

Rank: 2

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

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

本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

0

主题

3

帖子

67

积分

版主

Rank: 7Rank: 7Rank: 7

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

主题

40

帖子

206

积分

中级会员

Rank: 3Rank: 3

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

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

使用道具 举报

1

主题

16

帖子

160

积分

注册会员

Rank: 2

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

帖子

147

积分

注册会员

Rank: 2

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

主题

40

帖子

206

积分

中级会员

Rank: 3Rank: 3

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

帖子

147

积分

注册会员

Rank: 2

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

帖子

144

积分

注册会员

Rank: 2

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

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

使用道具 举报

0

主题

25

帖子

162

积分

注册会员

Rank: 2

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

主题

4

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 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-3-2 02:26 , Processed in 0.056302 second(s), 38 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.