搜索
楼主: Jimmy

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

  [复制链接]

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-1-9 23:44:03 | 显示全部楼层
本帖最后由 xxnku 于 2017-5-7 15:56 编辑

183-python-刘自军
按着视频的方法过了一遍:
[Python] 纯文本查看 复制代码
import re
import sys
from collections import OrderedDict
args = sys.argv
chr_dict = OrderedDict()

with open(args[1]) as hg19 :
        for line in hg19:
                line = line.rstrip()
                if line.startswith('>'):
                        temp_chr = line[1:]
                        chr_dict[temp_chr] = ''
                else:
                        chr_dict[temp_chr] += line.upper()
                        
for seqName,seq in chr_dict.items():
        #print (seqName + " =====> " + seq)
        seqLen = len(seq)
        N = seq.count("N")
        GC = seq.count("G") + seq.count("C") 
        print (seqName,seqLen, "%.2f"%(N/seqLen), "%.2f"%(GC/(seqLen-N)))


没看视频之前只会求一个chr的长度,憋了半天没憋出来如何求多个chr的长度,不知道如何存进dict里面,现在终于知道怎么把每个temp_name和line存到dict里面啦,哈哈哈。   而且这次视频录制效果提高了不少啊~~~~
最后关于视频里面的问题,如何是hg19的所有序列,该如何读取啊???用for line in hg19.readlines():   么???
回复 支持 反对

使用道具 举报

0

主题

5

帖子

36

积分

新手上路

Rank: 1

积分
36
发表于 2017-1-10 19:38:53 | 显示全部楼层
本帖最后由 死小孩 于 2017-1-14 09:22 编辑

109-R+perl+python
perl的思路是,先输出第一行,在匹配到“>”后,将上一个基因的结果输出,换行,输出这一个基因名。文件读完后,输出最后一个基因的数据。
[Perl] 纯文本查看 复制代码
#计算人类基因组每条染色体的长度和ATGC的总数
open UO,"$ARGV[0]" || die "Cannot open hg19_file $ARGV[0] \n";
open OUTFILE,">$ARGV[1]" || die "Cannot open output_file $ARGV[1] \n";

my $line = <UO>;
chomp $line;
print OUTFILE "$line"."\t";
while(<UO>){
        $line[$i] = $_;
        if ($line[$i] =~ /^\>/){
                chomp $line[$i];
                print OUTFILE "A:$A\t"."T:$T\t"."G:$G\t"."C:$C\t"."N:$N\t"."length:$chr_len\n"."$line[$i]\t";
                undef $A;    #清空$A;
                undef $T;
                undef $G;
                undef $C;
                undef $N;
                undef $chr_len;
        }
        else{
                $A += ($line[$i]=~tr/Aa//);
                $G += ($line[$i]=~tr/Gg//);
                $C += ($line[$i]=~tr/Cc//);
                $T += ($line[$i]=~tr/Tt//);
                $N += ($line[$i]=~tr/Nn//);
                $chr_len = $chr_len + $A + $T + $G + $C + $N;
        }
}
print OUTFILE "A:$A\t"."T:$T\t"."G:$G\t"."C:$C\t"."N:$N\t"."length:$chr_len\n";


Python

from collections import OrderedDict
chr_dict=OrderedDict()
temp_chr=''
seqlen=''
seq=''
outtxt=open('out.txt','w')
with open("hg19.fasta","r") as hg19:
    for line in hg19:
        line=line.strip()
        if line.startswith(">"):
            if seq!='':
                seqlen=len(seq)
                N=seq.count("N")
                GC=seq.count("G")+seq.count("g")+seq.count("c")+seq.count("C")
                print(temp_chr,seqlen,"%.2f"%(N/seqlen),"%.2f"%(GC/(seqlen-N)))
            temp_chr=line
            seq=''
        else:
            seq+=line
seqlen=len(seq)
N=seq.count("N")
GC=seq.count("G")+seq.count("g")+seq.count("c")+seq.count("C")
print(temp_chr,seqlen,"%.2f"%(N/seqlen),"%.2f"%(GC/(seqlen-N)))

回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-1-10 22:49:19 | 显示全部楼层
learnyoung 发表于 2017-1-7 15:03
下载了1号染色体的数据进行一个测试,我的电脑h19估计跑不动,但没想到的是chr1一个文件也跑了很久。
解题 ...

哥们,问你个问题哈。
python 中如果我想将print的结果不是直接输出到屏幕上,而是将print的结果输出到一个文件中,应该怎么办呢?
回复 支持 反对

使用道具 举报

1

主题

31

帖子

450

积分

中级会员

Rank: 3Rank: 3

积分
450
发表于 2017-1-10 23:59:32 | 显示全部楼层
本帖最后由 x2yline 于 2017-1-23 14:38 编辑

077 R-python-shell-x2yline
(1)R语言实现

数据来源:
demo数据 http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz

[AppleScript] 纯文本查看 复制代码
setwd('E:\\r\\biotrainee_demo\\class 2')
# 读入数据
remove(list=ls())
t1 <- Sys.time()
directory <-'chrY.fa' 
df <- read.csv(directory, header=F, stringsAsFactors=F)
# index_df 为chr所在的位置
index_df <- data.frame(begin=which(sapply(df[,1], function(x){
  substr(x, start=1, stop=1)=='>'})))
# index_df1 为string所在的位置+1
index_df1 <- data.frame(rbind(matrix(index_df[-1,1]),dim(df)[1]+1))
# 把index_start和index_end存入data.frame
index_df2 <- cbind(index_df, index_df1)
remove(index_df1, index_df)
# 得出每个染色体对应string后计算其N与GC百分比
result <- apply(index_df2, 1, function(x) {  # 把提取字符串后把字符串变为大写
  y <- toupper(paste(df[(x[1]+1):(x[2]-1),1], collapse=''))
  y <- matrix(strsplit(y, split=character(0))[[1]])
  N <- length(which(y[,1] =='N'))
  all_base_num <- dim(y)[1]
  N_ratio <- N/all_base_num 
  GC <- length(which(y[,1] =='G' | y[,1] == 'C'))
  GC_ratio <- GC/(all_base_num - N)
  c(N, N_ratio, GC, GC_ratio, all_base_num)
})
# 把行名改为N和GC并转秩
rownames(result) = c('N', 'N_ratio', 'GC', 'GC_ratio', 'all_base_num')
result <- t(result)
# 取结果前几行
head(result)
difftime(Sys.time(), t1, units = 'secs')


由于电脑问题,试了一下1号染色体,电脑卡住了,于是又试了一下Y染色体,跑出来结果如下:
耗时:Time difference of 41.44945 secs
            N                N_ratio         GC              GC_ratio         all_base_num
>chrY 33720000   0.5679295   10252459   0.3996504     59373566

电脑配置信息:
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

(2)python3实现
数据来源: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
数据下载时间:2017-01-10 23:08

运行结果:
chr        GC_ratio        N_ratio        Length        N        A        T        C        G
>chr1        0.41744        0.09617        249250621        23970000        65570891        65668756        47024412        47016562
>chr2        0.40244        0.02054        243199373        4994855        71102632        71239379        47915465        47947042
>chr3        0.39694        0.01629        198022430        3225295        58713343        58760485        38653197        38670110
>chr4        0.38248        0.01827        191154276        3492600        57932980        57952068        35885806        35890822
>chr5        0.39516        0.0178        180915260        3220000        53672554        53804137        35089383        35129186
>chr6        0.39611        0.02174        171115067        3720001        50554433        50533923        33143287        33163423
>chr7        0.40751        0.02378        159138663        3785000        45997757        46047257        31671670        31636979
>chr8        0.40176        0.02374        146364022        3475100        42767293        42715025        28703983        28702621
>chr9        0.41317        0.14921        141213431        21070000        35260078        35243882        24826212        24813259
>chr10        0.41585        0.03114        135534747        4220009        38330752        38376915        27308648        27298423
>chr11        0.41566        0.02872        135006516        3877000        38307244        38317436        27236798        27268038
>chr12        0.40812        0.02518        133851895        3370502        38604831        38624517        26634995        26617050
>chr13        0.38527        0.17001        115169878        19580000        29336945        29425459        18412698        18414776
>chr14        0.40887        0.17755        107349540        19060000        25992966        26197495        18027132        18071947
>chr15        0.42201        0.20322        102531392        20836626        23620876        23597921        17247582        17228387
>chr16        0.44789        0.12694        90354753        11470000        21724083        21828642        17630040        17701988
>chr17        0.4554        0.04187        81195210        3400000        21159933        21206981        17727956        17700340
>chr18        0.39785        0.0438        78077248        3420019        22465380        22489493        14838685        14863671
>chr19        0.4836        0.05615        59128983        3320000        14390632        14428951        13478255        13511145
>chr20        0.44126        0.05585        63025520        3520000        16523053        16725227        13107828        13149412
>chr21        0.40833        0.27059        48129895        13023253        10422924        10348785        7160212        7174721
>chr22        0.47988        0.31985        51304566        16410021        9094775        9054551        8375984        8369235
>chrX        0.39496        0.02686        155270560        4170000        45648952        45772424        29813353        29865831
>chrY        0.39965        0.56793        59373566        33720000        7667625        7733482        5099171        5153288


运行消耗时间:309 seconds

代码如下:
[Python] 纯文本查看 复制代码
import os
import time
begin = time.time()
os.chdir(r'F:\tmp\chromFa')
def count_n_and_gc(file):
    content = []
    chromsome = []
    g = 0; c = 0; n = 0; a = 0; t = 0
    with open(file) as f:
        raw_list = f.readlines()
        for i in raw_list:
            if not i.startswith('>'):
                i = i.upper()
                n +=  i.count('N')
                g += i.count('G')
                c += i.count('C')
                a += i.count('A')
                t += i.count('T')
            else:
                if  chromsome:
                    content.append((n ,a, t, c, g))
                    g = 0; c = 0; n = 0; a = 0; t = 0
                chromsome.append(i.strip())
        content.append(( n ,a, t, c, g))
    return (content,chromsome)
content = []
chromsome = []
for i in (list(range(1,23)) + ['X','Y']):
    file = 'chr'+ str(i) + '.fa'
    print('Start dealing with ' + file)
    m, n = count_n_and_gc(file)
    content += m
    chromsome += n
all_info = 'chr,GC_ratio,N_ratio,Length,N,A,T,C,G'
for i in range(len(chromsome)):
    data = '\n'+str(chromsome) +',' + "%.5f"%((content[-1]+content[-2])/sum(content[1:])) +','  + "%.5f" %(content[0]/(sum(content))) +','  +str((sum(content))) +','  +str((content[0])) + ','  +str(content[1])+',' +str(content[2])+','  +str(content[3])+','  +str(content[4]) 
    all_info += data
with open('hg19_analysis.csv','w') as f:
    f.write(all_info)
print('Time using:'+ str(time.time() - begin) + ' seconds\n')

(3)shell +python3

在windows下使用shell:
按照如下博客在windows下搭建shell使用环境
http://www.hi-linux.com/2016/05/ ... 5%9E%E5%99%A8Babun/

先使用shell脚本把所有chromFa.tar.gz 中的所有.fa文件合并为一个hg19.fa文件
脚本如下
[Shell] 纯文本查看 复制代码
tar zvfx chromFa.tar.gz    
cat *.fa > hg19.fa  
rm chr*.fa  
less hg19.fa 
按照老师的方法对python算法进行改良
改良后的代码如下:
[Python] 纯文本查看 复制代码
import os
import time
import re
import sys
from collections import OrderedDict
start = time.clock()
def count_fasta_atcgn(file_path, buffer_size=1024*1024):
    bases = ['N', 'A', 'T', 'C', 'G']
    ATCG_analysis = OrderedDict()
    with open(file_path, 'r') as f:
        line1 = f.readline().upper()
        chr_i = re.split('\s', line1)[0][1:]
        print(chr_i)
        ATCG_analysis[chr_i] = OrderedDict()
        for base in bases:
            ATCG_analysis[chr_i][base] = 0
        while True:
            chunk = f.read(buffer_size).upper()
            if '>' in chunk:
                chromsome = re.split('>',chunk)
                if chromsome[0]:
                    for base in bases:
                        ATCG_analysis[chr_i][base] += chromsome[0].count(base)
                for i in chromsome[1:]:
                    if i:
                        chr_i = re.split('\s', i[0:i.index('\n')])[0]
                        print(chr_i)
                        strings_i = i[i.index('\n'):]
                        ATCG_analysis[chr_i] = OrderedDict()
                        for base in bases:
                            ATCG_analysis[chr_i][base] = strings_i.count(base)
            else:
                for base in bases:
                    ATCG_analysis[chr_i][base] += chunk.count(base)
            if not chunk:
                break
    return ATCG_analysis

def write_atcg_to_csv(ATCG_analysis, file_path = '.'):
    file = os.path.join(file_path,'atcg_analysis.csv')
    csv_content = 'chromsome\tGC_content\tN_content\tLength\tN\tA\tT\tC\tG\n'
    for chr_id, atcg_count in ATCG_analysis.items():
        GC = atcg_count['G'] + atcg_count['C']
        N = atcg_count['N']
        Length = sum(atcg_count.values())
        GC_content = GC*1.0/(Length-N)
        N_content = N*1.0/Length
        csv_content += chr_id + '\t' + '%.4f'%GC_content + '\t' + '%.4f'%N_content + '\t' + str(Length) + '\t' + str(atcg_count['N']) +'\t' + str(atcg_count['A']) + '\t' + str(atcg_count['T']) + '\t' + str(atcg_count['C'])+'\t'+ str(atcg_count['G'])+ '\n'
    with open(file, 'w') as f:
        csv_file_content = re.sub('\t', ',', csv_content)
        f.write(csv_file_content)
    print(u'File have been saved in '+ file)
    return csv_content

    
if sys.argv:
    result = OrderedDict()
    for f in sys.argv:
        done = 0
        f= f.strip(''''"''')
        if f[-2:] == 'py' or not os.path.exists(f):
            continue
        print(f)
        
        try:
            done = 1
            result = OrderedDict(count_fasta_atcgn(file_path = f, buffer_size = 1024*2048), **result)
        except Exception as e:
            if f.startswith('-'):
                pass
            else:
                print(type(e))
    
    if done == 1:
        file = write_atcg_to_csv(result)
        print(file)
        print('used %.2f s'%(time.clock()-start))
    else:
        print ('\n\nSorry! The command is invalid!\n')
else:
    directory = input('Enter your file: ')
    start = time.clock()
    if directory.count('.') != 1 or directory[-2:] == 'py' or not os.path.exists(directory):
        print('Your file is invalid!')
    else:
        result = count_fasta_atcgn(file_path = directory, buffer_size = 1024*2048)
        file = write_atcg_to_csv(result)
        print('used %.2f s'%(time.clock()-start))

保存代码为fasta_atcgn_summary.py 文件


在命令行下输入:
python fasta_atcgn_summary.py   F:\tmp\hg19.fa

输出结果为:
chromsome        GC_content        N_content        Length        N        A        T        C        G
CHRUN_GL000237        0.4666        0        45867        0        12273        12191        10241        11162
CHR14        0.4089        0.1776        107349540        19060000        25992966        26197495        18027132        18071947
CHR9_GL000199_RANDOM        0.3791        0        169874        0        54702        50765        34981        29426
CHR2        0.4024        0.0205        243199373        4994855        71102632        71239379        47915465        47947042
CHR8_GL000197_RANDOM        0.5401        0.0027        37175        100        8644        8408        9883        10140
CHRUN_GL000247        0.436        0        36422        0        11002        9540        7794        8086
CHR19        0.4836        0.0561        59128983        3320000        14390632        14428951        13478255        13511145
CHR15        0.422        0.2032        102531392        20836626        23620876        23597921        17247582        17228387
CHR17_GL000206_RANDOM        0.5309        0        41001        0        8770        10463        10483        11285
CHRY        0.3997        0.5679        59373566        33720000        7667625        7733482        5099171        5153288
CHRUN_GL000226        0.3903        0        15008        0        4502        4649        2626        3231
CHRM        0.4449        0        16571        0        5113        4086        5192        2180
CHR7        0.4075        0.0238        159138663        3785000        45997757        46047257        31671670        31636979
CHRUN_GL000232        0.4182        0        40652        0        11534        12116        8490        8512
CHR9_GL000201_RANDOM        0.5944        0        36148        0        7101        7560        10373        11114
CHRX        0.395        0.0269        155270560        4170000        45648952        45772424        29813353        29865831
CHR9        0.4132        0.1492        141213431        21070000        35260078        35243882        24826212        24813259
CHRUN_GL000245        0.3631        0        36651        0        12875        10467        6707        6602
CHRUN_GL000227        0.4101        0        128374        0        41076        34652        26250        26396
CHRUN_GL000214        0.4152        0        137718        0        40645        39891        27484        29698
CHR6_SSTO_HAP7        0.4421        0.1532        4928567        755016        1172219        1156168        918220        926944
CHRUN_GL000239        0.454        0        33824        0        9221        9246        7268        8089
CHR4_CTG9_HAP1        0.3638        0        590426        0        185171        190487        107556        107212
CHR13        0.3853        0.17        115169878        19580000        29336945        29425459        18412698        18414776
CHR4_GL000194_RANDOM        0.4326        0        191469        0        52949        55693        41521        41306
CHRUN_GL000225        0.4765        0        211173        0        56943        53599        48931        51700
CHR1_GL000191_RANDOM        0.4435        0        106433        0        27971        31264        23785        23413
CHR17_GL000205_RANDOM        0.4173        0        174588        0        49255        52471        36716        36146
CHRUN_GL000230        0.4171        0        43691        0        12678        12789        9338        8886
CHR6_MANN_HAP4        0.4415        0.1244        4683263        582522        1156084        1134193        903022        907442
CHR4_GL000193_RANDOM        0.428        0        189789        0        53509        55056        40616        40608
CHRUN_GL000243        0.4601        0        43341        0        10728        12673        10304        9636
CHR3        0.3969        0.0163        198022430        3225295        58713343        58760485        38653197        38670110
CHR18_GL000207_RANDOM        0.4395        0        4262        0        1004        1385        588        1285
CHR7_GL000195_RANDOM        0.4066        0        182896        0        53469        55057        37021        37349
CHR20        0.4413        0.0559        63025520        3520000        16523053        16725227        13107828        13149412
CHR6_DBB_HAP3        0.4493        0.0881        4610396        406094        1172535        1142951        943550        945266
CHR11_GL000202_RANDOM        0.5461        0        40103        0        9226        8978        11254        10645
CHRUN_GL000219        0.3996        0        179198        0        54530        53059        35501        36108
CHR6_QBL_HAP6        0.4509        0.0687        4611984        316659        1189669        1168711        967043        969902
CHRUN_GL000231        0.4467        0        27386        0        7100        8054        6203        6029
CHRUN_GL000213        0.409        0        164239        0        48047        49015        33831        33346
CHRUN_GL000244        0.4363        0        39929        0        10948        11560        8564        8857
CHRUN_GL000238        0.4        0        39939        0        10404        13559        7805        8171
CHRUN_GL000228        0.5394        0        129120        0        30512        28967        35495        34146
CHRUN_GL000248        0.4566        0        39786        0        10114        11507        9265        8900
CHR9_GL000200_RANDOM        0.3995        0        187035        0        55353        56966        37202        37514
CHR1_GL000192_RANDOM        0.4149        0        547496        0        163078        157247        112730        114441
CHRUN_GL000218        0.4162        0        161147        0        46030        48047        33296        33774
CHR6_MCF_HAP5        0.4532        0.2149        4833398        1038487        1050968        1024227        858581        861135
CHR9_GL000198_RANDOM        0.3786        0        90085        0        27843        28140        15617        18485
CHRUN_GL000211        0.3871        0        166566        0        50926        51165        31968        32507
CHR8        0.4018        0.0237        146364022        3475100        42767293        42715025        28703983        28702621
CHR21_GL000210_RANDOM        0.543        0.0036        27682        100        6288        6317        7726        7251
CHR19_GL000209_RANDOM        0.465        0        159169        0        43978        41183        35443        38565
CHR6_COX_HAP2        0.4468        0        4795371        0        1341236        1311570        1069806        1072759
CHRUN_GL000234        0.4307        0        40531        0        10797        12277        8725        8732
CHR17        0.4554        0.0419        81195210        3400000        21159933        21206981        17727956        17700340
CHR10        0.4158        0.0311        135534747        4220009        38330752        38376915        27308648        27298423
CHR18        0.3978        0.0438        78077248        3420019        22465380        22489493        14838685        14863671
CHR5        0.3952        0.0178        180915260        3220000        53672554        53804137        35089383        35129186
CHRUN_GL000215        0.42        0        172545        0        50334        49738        36250        36223
CHRUN_GL000220        0.4846        0        161802        0        37230        46155        40720        37697
CHR19_GL000208_RANDOM        0.3766        0        92689        0        29179        28602        18369        16539
CHR6        0.3961        0.0217        171115067        3720001        50554433        50533923        33143287        33163423
CHRUN_GL000240        0.4255        0        41933        0        13601        10490        8962        8880
CHR17_CTG5_HAP1        0.4546        0.0595        1680828        100000        429214        432922        355909        362783
CHRUN_GL000235        0.38        0        34474        0        11845        9530        6585        6514
CHR17_GL000203_RANDOM        0.3329        0        37498        0        12564        12452        6074        6408
CHRUN_GL000249        0.4678        0        38502        0        10793        9698        8978        9033
CHR4        0.3825        0.0183        191154276        3492600        57932980        57952068        35885806        35890822
CHRUN_GL000223        0.432        0        180455        0        52931        49568        38849        39107
CHR11        0.4157        0.0287        135006516        3877000        38307244        38317436        27236798        27268038
CHRUN_GL000217        0.3759        0        172149        0        51945        55495        32578        32131
CHRUN_GL000229        0.5015        0        19913        0        3983        5944        5385        4601
CHRUN_GL000222        0.4397        0        186861        0        51533        53158        40866        41304
CHR8_GL000196_RANDOM        0.3965        0        38914        0        13843        9642        7910        7519
CHR16        0.4479        0.1269        90354753        11470000        21724083        21828642        17630040        17701988
CHR22        0.4799        0.3199        51304566        16410021        9094775        9054551        8375984        8369235
CHRUN_GL000233        0.4242        0        45941        0        12192        14260        9504        9985
CHR1        0.4174        0.0962        249250621        23970000        65570891        65668756        47024412        47016562
CHRUN_GL000216        0.4197        0        172294        0        41409        58566        46717        25602
CHRUN_GL000224        0.4329        0        179693        0        50248        51660        37430        40355
CHRUN_GL000242        0.484        0        43523        0        10591        11869        10033        11030
CHRUN_GL000221        0.3864        0        155397        0        47915        47444        29886        30152
CHR6_APD_HAP1        0.4393        0.4979        4622290        2301543        660350        640840        509167        510390
CHRUN_GL000236        0.4163        0        41934        0        13006        11470        8432        9026
CHRUN_GL000212        0.4438        0        186858        0        52729        51193        42454        40482
CHRUN_GL000241        0.3729        0        42152        0        13620        12815        7864        7853
CHR17_GL000204_RANDOM        0.5447        0        81310        0        19702        17322        22701        21585
CHR12        0.4081        0.0252        133851895        3370502        38604831        38624517        26634995        26617050
CHRUN_GL000246        0.3865        0        38154        0        10440        12968        7354        7392
CHR21        0.4083        0.2706        48129895        13023253        10422924        10348785        7160212        7174721


used 166.18 s

在命令行下输入:
python fasta_atcgn_summary.py   demo.txt

输出结果为:
chromsome       GC_content      N_content       Length  N       A       T       C       G
CHR_4   0.4091  0.1200  25      3       9       4       2       7
CHR_3   0.4483  0.1212  33      4       10      6       3       10
CHR_2   0.4333  0.1176  34      4       11      6       5       8
CHR_1   0.4250  0.0909  44      4       13      10      7       10

used 0.03 s

python可视化处理
[Python] 纯文本查看 复制代码
import os
import time
import re
import sys
from collections import OrderedDict
start = time.clock()
def count_fasta_atcgn(file_path, buffer_size=1024*1024):
    bases = ['N', 'A', 'T', 'C', 'G']
    ATCG_analysis = OrderedDict()
    with open(file_path, 'r') as f:
        line1 = f.readline().upper()
        chr_i = re.split('\s', line1)[0][1:]
        print(chr_i)
        ATCG_analysis[chr_i] = OrderedDict()
        for base in bases:
            ATCG_analysis[chr_i][base] = 0
        while True:
            chunk = f.read(buffer_size).upper()
            if '>' in chunk:
                chromsome = re.split('>',chunk)
                if chromsome[0]:
                    for base in bases:
                        ATCG_analysis[chr_i][base] += chromsome[0].count(base)
                for i in chromsome[1:]:
                    if i:
                        chr_i = re.split('\s', i[0:i.index('\n')])[0]
                        print(chr_i)
                        strings_i = i[i.index('\n'):]
                        ATCG_analysis[chr_i] = OrderedDict()
                        for base in bases:
                            ATCG_analysis[chr_i][base] = strings_i.count(base)
            else:
                for base in bases:
                    ATCG_analysis[chr_i][base] += chunk.count(base)
            if not chunk:
                break
    return ATCG_analysis
 
def write_atcg_to_csv(ATCG_analysis, file_path = '.'):
    file = os.path.join(file_path,'atcg_analysis.csv')
    csv_content = 'chromsome\tGC_content\tN_content\tLength\tN\tA\tT\tC\tG\n'
    for chr_id, atcg_count in ATCG_analysis.items():
        GC = atcg_count['G'] + atcg_count['C']
        N = atcg_count['N']
        Length = sum(atcg_count.values())
        GC_content = GC*1.0/(Length-N)
        N_content = N*1.0/Length
        csv_content += chr_id + '\t' + '%.4f'%GC_content + '\t' + '%.4f'%N_content + '\t' + str(Length) + '\t' + str(atcg_count['N']) +'\t' + str(atcg_count['A']) + '\t' + str(atcg_count['T']) + '\t' + str(atcg_count['C'])+'\t'+ str(atcg_count['G'])+ '\n'
    with open(file, 'w') as f:
        csv_file_content = re.sub('\t', ',', csv_content)
        f.write(csv_file_content)
    print(u'File have been saved in '+ file)
    return csv_content
    
file_path = 'F:\genome\chromFa\hg19.fa'


ATCG_analysis = count_fasta_atcgn(file_path, buffer_size=1024*1024)
cg_list = []
chr_id_list = list(range(1,23)) + ['X','Y','M']
for i in chr_id_list:
    cg_list.append((ATCG_analysis['CHR'+str(i)]['G']+ATCG_analysis['CHR'+str(i)]['C'])/(ATCG_analysis['CHR'+str(i)]['A']+ATCG_analysis['CHR'+str(i)]['T']+ATCG_analysis['CHR'+str(i)]['C']+ATCG_analysis['CHR'+str(i)]['G'])*100)
import matplotlib.pyplot as plt
plt.bar(left = range(25), height = cg_list, color='k')
for i in range(len(cg_list)):
    plt.text( x=i -.1, y=cg_list[i]+.35,s=str(round(cg_list[i])))
plt.title('GC content for hg19 genome')
plt.ylabel('GC content (%)')
pos = []
for i in range(len(chr_id_list)):
    pos.append(i + 0.35)
plt.xticks(pos, list(range(1,23)) + ['X','Y','MT'], fontsize=8)
plt.xlim(-0.2, )
plt.ylim(0, 100)
plt.savefig('F:\gc.png',dpi=600)
plt.show()






回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

84

积分

注册会员

Rank: 2

积分
84
发表于 2017-1-11 10:04:17 | 显示全部楼层
042-r+shell+python-鹏程  按嘉宾视频中的来,跟着一点点模仿学习,

[Python] 纯文本查看 复制代码
# coding: utf-8

# In[3]:

import sys 
import os


# In[59]:

argvs=sys.argv
d={}
bases = ("A","T","G","C","N")
#with open("testdata.txt","rt") as f:
with open(argvs[1],"rt") as f:
    for line in f:
        line=line.strip()##先把空格换行去掉
        if line.startswith(">"):
            chr_id=line[1:]
            d[chr_id]={}
            for base in bases:
                d[chr_id][base]=0
        else:
            line=line.upper()
            for base in bases:
                d[chr_id][base]+=line.count(base)
for chr_id, value in d.items():
    print (chr_id)
    length=sum(d[chr_id].values())
    GC=d[chr_id]["G"]+d[chr_id]["C"]
    N=d[chr_id]["N"]
    #print(length,GC)
    print ("chromesome length is: %d" % length)
    print ("N content is: %.2f"%(N/length))
    print ("GC content is: %.2f"%(GC/length))
    #print ("T:",d[chr_id]["T"])
    #print ("ALL:",d[chr_id]["A"]/d[chr_id]["T"] )
    #print (value)
    
          


# In[ ]:



回复 支持 反对

使用道具 举报

633

主题

1172

帖子

3947

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3947
 楼主| 发表于 2017-1-11 14:37:24 | 显示全部楼层
x2yline 发表于 2017-1-10 23:59
077 R-python-shell-x2yline
(1)R语言实现

如果你想用R语言解决这个问题,请试一试biostring这个包
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

1

主题

39

帖子

355

积分

中级会员

Rank: 3Rank: 3

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

使用道具 举报

1

主题

39

帖子

355

积分

中级会员

Rank: 3Rank: 3

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

我试用了下这个包 很好用。。
人生若只如初见!
回复 支持 反对

使用道具 举报

1

主题

41

帖子

270

积分

中级会员

Rank: 3Rank: 3

积分
270
发表于 2017-1-11 17:07:58 | 显示全部楼层
xxnku 发表于 2017-1-10 22:49
哥们,问你个问题哈。
python 中如果我想将print的结果不是直接输出到屏幕上,而是将print的结果输出到一 ...

直接write进去一个文本文件就好了比如
[Python] 纯文本查看 复制代码
g=open('test.txt','w')
a=1
g.write(a)
g.close()
回复 支持 反对

使用道具 举报

1

主题

41

帖子

270

积分

中级会员

Rank: 3Rank: 3

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

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

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-12-19 06:11 , Processed in 0.116080 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.