搜索
查看: 1386|回复: 102

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

  [复制链接]

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

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

主题

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

使用道具 举报

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()
回复 支持 2 反对 0

使用道具 举报

7

主题

25

帖子

163

积分

注册会员

Rank: 2

积分
163
发表于 7 天前 | 显示全部楼层
本帖最后由 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

主题

4

帖子

55

积分

注册会员

Rank: 2

积分
55
发表于 7 天前 | 显示全部楼层
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

帖子

38

积分

新手上路

Rank: 1

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

主题

9

帖子

102

积分

注册会员

Rank: 2

积分
102
发表于 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-22 03:31 , Processed in 0.217973 second(s), 38 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.