搜索
查看: 5581|回复: 138

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

  [复制链接]

608

主题

1073

帖子

3577

积分

管理员

Rank: 9Rank: 9Rank: 9

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

主题

30

帖子

333

积分

中级会员

Rank: 3Rank: 3

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

使用道具 举报

8

主题

40

帖子

340

积分

版主

Rank: 7Rank: 7Rank: 7

积分
340
QQ
发表于 2017-1-15 17:31:15 | 显示全部楼层
本帖最后由 旭日早升 于 2017-3-11 21:44 编辑

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










回复 支持 4 反对 0

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

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

使用道具 举报

1

主题

34

帖子

253

积分

中级会员

Rank: 3Rank: 3

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

主题

7

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 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
回复 支持 2 反对 0

使用道具 举报

0

主题

3

帖子

63

积分

注册会员

Rank: 2

积分
63
发表于 2017-2-14 15:29:02 | 显示全部楼层
ID:179
Python
第二题确实学到了不少,思路一就是运用字典的方法来进行存储数据,然后提取字典的键和值来进行分析,但会很耗时间。另外一个就是用pysam包。为了解决第一题,电脑安装了64位python3版本的anaconda,使用conda list命令后,发现并没有pysam包。之后使用pip install pysam的命令安装pysam,结果也出现了安装到python2.7的路径下。用vim ~/.bashrc打开文件,在python2.7的路径前加上#后保存退出,再使用source ~/.bashrc进行重置。最后再使用pip install pysam就成功了。最后pysam安装在anaconda3/bin/目录下。pysam确实很好用。
[Python] 纯文本查看 复制代码
import pysam

hg19 = pysam.Fastafile('hg19.fa')
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,seqlength,'%.2f'%(N/seqlength),"%.2f"%(GC/(seqlength-N)))
回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

67

积分

注册会员

Rank: 2

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

使用道具 举报

14

主题

67

帖子

491

积分

中级会员

Rank: 3Rank: 3

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

主题

4

帖子

68

积分

注册会员

Rank: 2

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

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

本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

1

主题

16

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 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
回复 支持 反对

使用道具 举报

7

主题

34

帖子

192

积分

注册会员

Rank: 2

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

主题

41

帖子

250

积分

中级会员

Rank: 3Rank: 3

积分
250
发表于 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'))
回复 支持 反对

使用道具 举报

7

主题

34

帖子

192

积分

注册会员

Rank: 2

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

已更新代码与结果
欢迎关注博客: 生信客部落
回复 支持 反对

使用道具 举报

0

主题

12

帖子

101

积分

注册会员

Rank: 2

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

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

11

主题

44

帖子

210

积分

中级会员

Rank: 3Rank: 3

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

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

使用道具 举报

1

主题

34

帖子

253

积分

中级会员

Rank: 3Rank: 3

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

主题

12

帖子

101

积分

注册会员

Rank: 2

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

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

使用道具 举报

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

本版积分规则

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

GMT+8, 2017-6-29 18:39 , Processed in 0.063058 second(s), 40 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.