搜索
楼主: Jimmy

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

  [复制链接]

1

主题

15

帖子

154

积分

注册会员

Rank: 2

积分
154
发表于 2017-2-3 17:50:38 | 显示全部楼层
本帖最后由 nickshi 于 2017-2-3 17:52 编辑

198 python +R
学习了东老师的视频,收益颇多,对字典的应用提升了

[mw_shl_code=python,true]from collections import OrderedDict
import time
sum_atcg = OrderedDict()
bases = ["A","T","C","G","N"]
sum_all = 0
start =time.clock()
with open('h19.txt','r') as f:
    for line in f:
        if line.startswith('>'):
            chr_id = line[1:]
            sum_atcg[chr_id] = {}
            for base in bases:
                sum_atcg[chr_id][base] = 0
        else:
            line = line.upper().strip()
            sum_all += len(line)   #加了一个计算了全部长度
            for base in bases:
                sum_atcg[chr_id][base] += line.count(base)
    for chr_id ,atcg_count in sum_atcg.items():
        sum_chr = sum(atcg_count.values())
        GC =round( ( atcg_count['G'] + atcg_count["C"])*1.0/sum_chr,2) #不加1.0发现计算结果等于0?还好在视频里学到了,要不然自己遇到要百度很久
        print GC
        print chr_id
        for base in bases:
            print '%s : %s'%(base,atcg_count[base])
        print 'SUM_chr : %s'%sum_chr
        print 'GC: %s'%GC
    print 'sum_all:%s'%sum_all
end = time.clock()
print 'used %s s'% (end-start)[/mw_shl_code]
回复 支持 反对

使用道具 举报

0

主题

2

帖子

47

积分

新手上路

Rank: 1

积分
47
发表于 2017-2-4 16:56:38 | 显示全部楼层
070-python+R+shell-panda

perl版本:
[mw_shl_code=perl,true]#!/usr/bin/perl -w
use strict;
@ARGV==1 || die "Usage: perl $0 <ref.fa>\n";
open FA,"<$ARGV[0]";
my (%per,$id);
while (<FA>){
        chop;
        if (/^>(\S+)/){
                $id=$1;
        }else{
                $per{$id}{'chr'}+=length($_);
                my @_gc=();
                my @_n=();
                push @_gc,$_=~/[G,C]/gi;
                push @_n,$_=~/N/gi;
                $per{$id}{'gc'}+=scalar (@_gc) ;
                $per{$id}{'n'}+=scalar (@_n);
        }
}
close FA;
foreach(sort keys %per){
        print "$_\t$per{$_}{'chr'}\t".(sprintf("%.4f",$per{$_}{'gc'}/$per{$_}{'chr'}))."\t".(sprintf("%.4f",$per{$_}{'n'}/$per{$_}{'chr'}))."\n";
}[/mw_shl_code]

结果
chr1        249250621        0.3773        0.0962
chr10        135534747        0.4029        0.0311
chr11        135006516        0.4037        0.0287
chr11_gl000202_random        40103        0.5461        0.0000
chr12        133851895        0.3978        0.0252
chr13        115169878        0.3198        0.1700
chr14        107349540        0.3363        0.1776
chr15        102531392        0.3362        0.2032
chr16        90354753        0.3910        0.1269
chr17        81195210        0.4363        0.0419
chr17_ctg5_hap1        1680828        0.4276        0.0595
chr17_gl000203_random        37498        0.3329        0.0000
chr17_gl000204_random        81310        0.5447        0.0000
chr17_gl000205_random        174588        0.4173        0.0000
chr17_gl000206_random        41001        0.5309        0.0000
chr18        78077248        0.3804        0.0438
chr18_gl000207_random        4262        0.4395        0.0000
chr19        59128983        0.4564        0.0561
chr19_gl000208_random        92689        0.3766        0.0000
chr19_gl000209_random        159169        0.4650        0.0000
chr1_gl000191_random        106433        0.4435        0.0000
chr1_gl000192_random        547496        0.4149        0.0000
chr2        243199373        0.3942        0.0205
chr20        63025520        0.4166        0.0559
chr21        48129895        0.2978        0.2706
chr21_gl000210_random        27682        0.5410        0.0036
chr22        51304566        0.3264        0.3199
chr3        198022430        0.3905        0.0163
chr4        191154276        0.3755        0.0183
chr4_ctg9_hap1        590426        0.3638        0.0000
chr4_gl000193_random        189789        0.4280        0.0000
chr4_gl000194_random        191469        0.4326        0.0000
chr5        180915260        0.3881        0.0178
chr6        171115067        0.3875        0.0217
chr6_apd_hap1        4622290        0.2206        0.4979
chr6_cox_hap2        4795371        0.4468        0.0000
chr6_dbb_hap3        4610396        0.4097        0.0881
chr6_mann_hap4        4683263        0.3866        0.1244
chr6_mcf_hap5        4833398        0.3558        0.2149
chr6_qbl_hap6        4611984        0.4200        0.0687
chr6_ssto_hap7        4928567        0.3744        0.1532
chr7        159138663        0.3978        0.0238
chr7_gl000195_random        182896        0.4066        0.0000
chr8        146364022        0.3922        0.0237
chr8_gl000196_random        38914        0.3965        0.0000
chr8_gl000197_random        37175        0.5386        0.0027
chr9        141213431        0.3515        0.1492
chr9_gl000198_random        90085        0.3786        0.0000
chr9_gl000199_random        169874        0.3791        0.0000
chr9_gl000200_random        187035        0.3995        0.0000
chr9_gl000201_random        36148        0.5944        0.0000
chrM        16571        0.4449        0.0000
chrUn_gl000211        166566        0.3871        0.0000
chrUn_gl000212        186858        0.4438        0.0000
chrUn_gl000213        164239        0.4090        0.0000
chrUn_gl000214        137718        0.4152        0.0000
chrUn_gl000215        172545        0.4200        0.0000
chrUn_gl000216        172294        0.4197        0.0000
chrUn_gl000217        172149        0.3759        0.0000
chrUn_gl000218        161147        0.4162        0.0000
chrUn_gl000219        179198        0.3996        0.0000
chrUn_gl000220        161802        0.4846        0.0000
chrUn_gl000221        155397        0.3864        0.0000
chrUn_gl000222        186861        0.4397        0.0000
chrUn_gl000223        180455        0.4320        0.0000
chrUn_gl000224        179693        0.4329        0.0000
chrUn_gl000225        211173        0.4765        0.0000
chrUn_gl000226        15008        0.3903        0.0000
chrUn_gl000227        128374        0.4101        0.0000
chrUn_gl000228        129120        0.5394        0.0000
chrUn_gl000229        19913        0.5015        0.0000
chrUn_gl000230        43691        0.4171        0.0000
chrUn_gl000231        27386        0.4467        0.0000
chrUn_gl000232        40652        0.4182        0.0000
chrUn_gl000233        45941        0.4242        0.0000
chrUn_gl000234        40531        0.4307        0.0000
chrUn_gl000235        34474        0.3800        0.0000
chrUn_gl000236        41934        0.4163        0.0000
chrUn_gl000237        45867        0.4666        0.0000
chrUn_gl000238        39939        0.4000        0.0000
chrUn_gl000239        33824        0.4540        0.0000
chrUn_gl000240        41933        0.4255        0.0000
chrUn_gl000241        42152        0.3729        0.0000
chrUn_gl000242        43523        0.4840        0.0000
chrUn_gl000243        43341        0.4601        0.0000
chrUn_gl000244        39929        0.4363        0.0000
chrUn_gl000245        36651        0.3631        0.0000
chrUn_gl000246        38154        0.3865        0.0000
chrUn_gl000247        36422        0.4360        0.0000
chrUn_gl000248        39786        0.4566        0.0000
chrUn_gl000249        38502        0.4678        0.0000
chrX        155270560        0.3844        0.0269
chrY        59373566        0.1727        0.5679

real        14m30.784s
user        14m13.534s
sys        0m16.811s

python版本

[mw_shl_code=python,true]#!/usr/bin/env python
# -*- coding:utf-8 -*-
from __future__ import division
from collections import defaultdict
percent=defaultdict(dict)

with open ('ref.fa','rt') as f:
        for line in f:
                line=line.rstrip()
                if line.startswith('>'):
                        id=line.strip('>')
                        percent[id]['len']=0
                        percent[id]['gc']=0
                        percent[id]['n']=0
                else:
                        percent[id]['len']+=len(line)
                        percent[id]['gc']+=line.count("G")
                        percent[id]['gc']+=line.count("g")
                        percent[id]['gc']+=line.count("C")
                        percent[id]['gc']+=line.count("c")
                        percent[id]['n']+=line.count("N")
                        percent[id]['n']+=line.count("n")
for m in percent:
        gc=percent[m]['gc']/percent[m]['len']
        n=percent[m]['n']/percent[m]['len']
        print ('%s\t%d\t%.4f\t%.4f'%(m,percent[m]['len'],gc,n))[/mw_shl_code]

python比较差,写的比较烂,测试数据的结果:

chr_1        44        0.3864        0.0909
chr_3        33        0.3939        0.1212
chr_2        34        0.3824        0.1176
chr_4        25        0.3600        0.1200
回复 支持 反对

使用道具 举报

0

主题

2

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-2-5 14:22:17 | 显示全部楼层
本帖最后由 M_____X_ 于 2017-2-5 14:53 编辑

054-python
主要学习和巩固了一些python的基础知识和一些细节问题,虽然花了一些时间,但是感觉很有收获。
[mw_shl_code=python,true]
# 方法一

from collections import OrderedDict

chir_dirt = OrderedDict()
temp_chr = ""

with open("hg19.fasta", "r") as hg19:
    for line in hg19:
        line = line.strip()  # str.strip()文本处理函数,用来处理文本行末空格和换行符
        if line.startswith(">"):  # str.startswith()文本处理函数,用来检查字符串是否以指定字串开头
            temp_chr = line              # 存储键
            chir_dirt[temp_chr] = ""  # dir[ ] =    向字典中添加值
        else:
            chir_dirt[temp_chr] += line  # += 向字典中添加,=会覆盖

for Sname, Sseq in chir_dirt.items():  # dirt.items()返回可遍历的(键, 值) 元组数组,并分别存入变量
    # print (Sname, Sseq)
    seqLen = len(Sseq)
    Nnum = Sseq.count("N")
    GC = Sseq.count("G") + Sseq.count("g") + Sseq.count("c") + Sseq.count("C")
    print(Sname, seqLen, "%.2f" %
          (Nnum / seqLen), "%.2f" % (GC / (seqLen - Nnum)))

# 方法2

import pysam  # 利用pysam模块读取fasta格式文件的功能

hg19 = pysam.FastaFile("hg19.fasta")

for SName in hg19.references:  # .references chromosome 名称
    Sseq = hg19.fetch(SName)  # 依据references名称获取1行或多行
    seqLen = len(Sseq)
    Nnum = Sseq.count("N")
    GC = Sseq.count("G") + Sseq.count("g") + \
        Sseq.count("c") + Sseq.count("C")
    print(SName, seqLen, "%.2f" %
          (Nnum / seqLen), "%.2f" % (GC / (seqLen - Nnum)))[/mw_shl_code]

回复 支持 反对

使用道具 举报

4

主题

50

帖子

327

积分

中级会员

Rank: 3Rank: 3

积分
327
发表于 2017-2-7 14:09:21 | 显示全部楼层
129-python/r-纸雪

python版本的讲师代码的步骤简单明了,先读取fasta文件存到内存,后对内存的每个纪录统计输出,讲师还介绍了pysam,之前没接触过,我用测试数据试了一下,结果如下:
python代码
[mw_shl_code=python,true]from collections import OrderedDict
chr_dict = OrderedDict()
temp_chr = ""

with open("test.fasta","r") as hg19:
    for line in hg19:
        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():
    seqlen = len(seq)
    N = seq.count('N')
    GC =seq.count('G')+seq.count('g')+seq.count('C')+seq.count('c')
    print(seqname,seqlen,'%.2f'%(N/seqlen),"%.2f"%(GC/(seqlen-N)))

'''
#pip install pysam
import pysam
hg = pysam.Fastafile("hg19.fasta")
#list(zip(hg19.references,hg19.lengths))
#hg19.fetch('chrM')
for seqname in hg19.references:
    seq = hg19.fetch(seqname)
    seqlength = len(seq)
    N = seq.count('N')
    GC =seq.count('G')+seq.count('g')+seq.count('C')+seq.count('c')
    print(seqname,seqlen,'%.2f'%(N/seqlen),"%.2f"%(GC/(seqlen-N)))

'''
#>chr_1 44 0.09 0.42
#>chr_2 34 0.12 0.43
#>chr_3 33 0.12 0.45
#>chr_4 25 0.12 0.41[/mw_shl_code]


R的代码只用测试数据试了一下,出了结果,代码如下:
R代码
[mw_shl_code=applescript,true]seq<-readLines("test.fasta")
is.anno<-regexpr("^>",seq,perl=T)
seq.anno<-seq[which(is.anno==1)]#chom name
seq.context<-seq[which(is.anno==-1)]#sequences

#记录每条序列内容的行数,便于后来分组
start<-which(is.anno==1) #chom 开始行号
end<-start[2:length(start)]-1 #chom 最后一行的行号
end<-c(end,length(seq))
distance<-end-start #chom 总共的行数
index<-1:length(start)
index<-rep(index,distance)
seqs<-tapply(seq.context, index, paste,collapse='')
seq.context<-as.character(seqs)
seq.len<-nchar(seq.context) #DNA序列长度

#计算ATCGN含量
a_count<-''
t_count<-''
c_count<-''
g_count<-''
n_count<-''
gc_count<-''
gc_percent<-''
n_percent<-''

for(i in 1:length(seq.context))
{
  a_count<-length(gregexpr("[Aa]",seq.context,perl=T)[])
  t_count<-length(gregexpr("[Tt]",seq.context,perl=T)[])
  c_count<-length(gregexpr("[Cc]",seq.context,perl=T)[])
  g_count<-length(gregexpr("[Gg]",seq.context,perl=T)[])
  n_count<-length(gregexpr("[Nn]",seq.context,perl=T)[])
  gc_percent<-(as.numeric(g_count)+as.numeric(c_count))/as.numeric(seq.len)
  n_percent<-as.numeric(n_count)/as.numeric(seq.len)
}
result<-data.frame(seq.anno,seq.len,a_count,t_count,c_count,g_count,n_count,gc_percent,n_percent)#结果输出
print(result)
#  seq.anno seq.len a_count t_count c_count g_count n_count        gc_percent          n_percent
#1   >chr_1      44      13      10       7      10       4 0.386363636363636 0.0909090909090909
#2   >chr_2      34      11       6       5       8       4 0.382352941176471  0.117647058823529
#3   >chr_3      33      10       6       3      10       4 0.393939393939394  0.121212121212121
#4   >chr_4      25       9       4       2       7       3              0.36               0.12[/mw_shl_code]
思考:考虑到hg19基因组有几个G,若计算机内存很小的话,用python/r做的话,可以借用迭代器的思想,即读取一条完整纪录后,直接进行统计,以下一条纪录的“>”作为纪录之间的分割(最后一条记录除外),hg19的数据还在下,下完后用服务器(学生用的,只有1G内存)试试看。


回复 支持 反对

使用道具 举报

0

主题

5

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-2-7 23:01:45 | 显示全部楼层
本帖最后由 chzhniu 于 2017-2-7 23:05 编辑

182-perl
以测试文件进行脚本校正和测试。
1.png
思路:对测试文件进行逐行读取,然后对每一次读取进行判断,如果读取的是染色体“>XXX”,则输出染色体名称,否则分别进行ATGC含量统计。最后将统计结果输出。
初次写的脚本如下:
while(<>){
#chomp;#此处chomp之后结果不正确,可能是复制测试文件有误
if(/^>/){
       print;
       } else {
              chomp;
              $A =$_ =~tr/A|a/A/d;   #tr/Aa//,tr/Aa/A/效果一样,/d表示将匹配的进行替换
              $T =$_ =~tr/T|t/T/d;
              $G =$_ =~tr/G|g/G/d;
              $C =$_ =~tr/C|c/C/d;
              print"$A\t$T\t$G\t$C\t";
              }
}
输出结果为:
2.png
这个结果是对每一行进行统计然后将结果输出,并没有统计一条染色体中每种碱基总共含量。所以还需要一个累加过程;对脚本进行修改:
$A+=$_ =~tr/A|a/A/d;
结果如下:
3.png
输出结果为前边碱基的累加和。(file:///C:/Users/ADMINI~1/AppData/Local/Temp/msohtmlclip1/01/clip_image006.png试了很多方法就是搞不好。。。。。。)
以下是根据讲课内容修改的脚本:

[mw_shl_code=perl,true]while (<>){
chomp;
if (/^>/){
        $chr=$_;
        #print "$chr\n"
        }
        else {
                #chomp;
                $A{$chr} +=($_ =~tr/Aa//);
                $T{$chr} +=($_ =~tr/Tt//);
                $G{$chr} +=($_ =~tr/Gg//);
                $C{$chr} +=($_ =~tr/Cc//);
                $N{$chr} +=($_ =~tr/Nn//);
                }
}
foreach (sort keys %A){  #此处很关键,目的是取出染色体名称(键),后续根据键值计算碱基数量
$length=$A{$_}+$T{$_}+$G{$_}+$C{$_}+$N{$_};
$GC=($G{$_}+$C{$_})/$length;
print "$_\n$A{$_}\t$T{$_}\t$G{$_}\t$C{$_}\t$N{$_}\t$length\t$GC\n";
}[/mw_shl_code]

运行结果如下:
4.png
总结:本题目的关键是将文件转化成哈希,哈希键为染色体名称,值为每种碱基的数量。最后统计的时候,计算每一个键对应的各个哈希的值。同时,为是结果按染色体顺序排列,对哈希键进行sort。通过本体学会了碱基数量统计方法,同时学会了对哈希的转化。




回复 支持 反对

使用道具 举报

0

主题

2

帖子

18

积分

新手上路

Rank: 1

积分
18
发表于 2017-2-8 17:31:15 | 显示全部楼层
本帖最后由 颤抖的心丶 于 2017-2-8 17:33 编辑

195-perl
QQ截图20170208172714.png
这是结果:
QQ截图20170208172912.png
回复 支持 反对

使用道具 举报

0

主题

2

帖子

16

积分

新手上路

Rank: 1

积分
16
发表于 2017-2-9 16:53:02 | 显示全部楼层
电脑配置不行呀,拿出了老师视频讲解的几行数据进行了练习。学习了tapply,rep等函数的用法。理解了生物学知识,离自己上手还有一段距离。
[mw_shl_code=applescript,true]
setwd("E:/R")
seq<-readLines("hg19_demo.txt") # 读入序列,每个元素存入一行
is.anno<-regexpr("^>",seq,perl = T) # 正则匹配(regular expression)注释行,是注释行为1,否则为-1,然后根据1,-1来定义是染色体还是dna序列
mychr<-seq[which(is.anno==1)] # 选出染色体
myseq<-seq[which(is.anno==-1)] # 选出DNA系列
start<-which(is.anno==1) # 注释每条染色体所在的行数
end<-start[2:length(start)]-1  # 第二条记录注释行到最后一条记录注释行行号减一,即为每条记录结束行号,这里会统计少一行——最后一行的结束未统计
end <- c(end, length(seq) ) # 末尾添加一行:所有序列结束行
index <- 1:length(start) #start的长度 4个,分成4组
index <- rep(index, distance) # 分组标签
seqs <- tapply(seq.content, index, paste, collapse="") # 按照分组整合,拼接每条序列内容,返回一个列表,列表每个元素为一条序列的内容
seq.content<-as.character( seqs ) # 将列表转换为向量,向量每个元素为一条序列的内容
seq.len <- nchar(seq.content) # 获得序列长度
#计算ATCGN含量#
a_count<-""
t_count<-""
c_count<-""
g_count<-""
n_count<-""
gc_percent<-""
n_percent<-""
#for循环,正则匹配,
for(i in 1:length(seq.content)){
  a_count<-length(gregexpr("[Aa]",seq.content,perl=T)[])
  t_count<-length(gregexpr("[Tt]",seq.content,perl=T)[])
  c_count<-length(gregexpr("[Cc]",seq.content,perl=T)[])
  g_count<-length(gregexpr("[Gg]",seq.content,perl=T)[])
  n_count<-length(gregexpr("[Nn]",seq.content,perl=T)[])
  gc_percent<-(as.numeric(g_count)+as.numeric(c_count))/as.numeric(seq.len)
  n_percent<-as.numeric(n_count)/as.numeric(seq.len)
}
result<-data.frame(mychr,seq.len,a_count,t_count,c_count,g_count,n_count,gc_percent,n_percent)
result[/mw_shl_code]
回复 支持 反对

使用道具 举报

0

主题

15

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-2-9 17:32:23 | 显示全部楼层
本帖最后由 Aiyawq 于 2017-2-9 17:38 编辑

207-perl-wangQ;
这道题先拿测试数据试了一下,能够很快跑出结果
[mw_shl_code=perl,true]#!/usr/bin/perl
use strict;
use warnings;
open A,$ARGV[0] or die;
#open B,">",$ARGV[1] or die;
my ($chr,$seq);
my %hash;
my $time = `date`;
print "chrome\tcount_A\tcount_T\tcount_C\tcount_G\tcount_N\tlength\tpercent_GC\n";
while(<A>){
      chomp;
      if ($_ =~ /^>(.*)/){
      $chr = $1;
      $seq = "";
      }else{
            $seq .= $_;
            $hash{$chr} = $seq;
      }
}
close A;
foreach (sort keys %hash){
          my $len = length($hash{$_});
          my $count_A = $hash{$_} =~ tr/Aa//;
          my $count_T = $hash{$_} =~ tr/Tt//;
          my $count_C = $hash{$_} =~ tr/Cc//;
          my $count_G = $hash{$_} =~ tr/Gg//;
          my $count_N = $hash{$_} =~ tr/Nn//;
          my $percent_GC = ($count_C+$count_G)*100/($count_A+$count_T+$count_C+$count_G);
          print B "$_\t$count_A\t$count_T\t$count_C\t$count_G\t$count_N\t$len\t$percent_GC%\n";
}
print "$time\n";
print system ("wc -l $ARGV[0]");
#close B;[/mw_shl_code]
由于这个方法其实的效率并不高,因此跑其他染色体一时半会并没有跑出来,所以我只截取了测试数据和线粒体DNA数据的;后期技能提升了的话我想写一个效率高一点的重新跑出来分享给大家



这是测试数据的

这是测试数据的

这是线粒体的

这是线粒体的
回复 支持 反对

使用道具 举报

0

主题

3

帖子

47

积分

新手上路

Rank: 1

积分
47
发表于 2017-2-10 10:08:50 | 显示全部楼层
057-perl+shell-Sylvia
原先账号密码忘记了,新注册了个账号。基础太差,拖了好久,又是各种参考才做出来, hg19基因组按群主的步骤下载了。
[mw_shl_code=perl,true]#!/usr/bin/perl
my $hg19 = $ARGV[0];
open IN,$hg19 or die "$!";
$/=">",<IN>;
print "chr\tnumber_a\tnumber_t\tnumber_c\tnumber_g\tnumber_n\tLen\tGC(%)\tN(%)\n";
while(<IN>){
chomp;
my($chr,$seq)=split/\n/,$_,2;
  foreach ($seq){
  chomp;
  $a = ($_=~tr/Aa//);
  $c = ($_=~tr/Cc//);
  $g = ($_=~tr/Gg//);
  $t = ($_=~tr/Tt//);
  $n = ($_=~tr/Nn//);
  $length = $a+$c+$g+$t+$n;
  $GC = (($g+$c)/($a+$c+$g+$t))*100;
  $N = ($n/$length)*100;

  printf"%s\t%d\t%d\t%d\t%d\t%d\t%d\t%.2f\t%.2f\n",$chr,$a,$t,$c,$g,$n,$length,$GC,$N;
}
}
[/mw_shl_code]
回复 支持 反对

使用道具 举报

0

主题

3

帖子

53

积分

注册会员

Rank: 2

积分
53
发表于 2017-2-10 20:33:36 | 显示全部楼层
本帖最后由 Aliangkou 于 2017-2-10 20:35 编辑

145-R-翠湖心影


[mw_shl_code=python,true]import os
from collections import OrderedDict

chr_dict = OrderedDict()
temp_chr = ""

os.chdir('/Users/Administrator/Desktop/python') #工作目录
with open("hg19.fasta","rt") as f:
        for line in f:
                line = line.strip()
                if line.startswith(">"):
                        for Name,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")
                                GC = C + G
                                sum = A + T + C + G + N
                                print(Name)
                                print("A: %s" % (A))
                                print("G: %s" % (G))
                                print("C: %s" % (C))
                                print("T: %s" % (T))
                                print("N: %s" % (N))
                                print("SUM:  %s" % (sum))
                                print("GC: %s" % (GC/sum))
                        temp_chr = line
                        chr_dict = OrderedDict()
                        chr_dict[temp_chr] = ""
                else:
                        chr_dict[temp_chr] += line

for Name,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")
        GC = C + G
        sum = A + T + C + G + N
        print(Name)
        print("A: %s" % (A))
        print("G: %s" % (G))
        print("C: %s" % (C))
        print("T: %s" % (T))
        print("N: %s" % (N))
        print("SUM:  %s" % (sum))
        print("GC: %s" % (GC/sum))
       [/mw_shl_code]


C:\Users\Administrator\Desktop\python\20170210203209.png
20170210203209.png
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-2-29 07:48 , Processed in 0.031206 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.