搜索
查看: 31574|回复: 143

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

  [复制链接]

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 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文件的这些探究。
[mw_shl_code=applescript,true]>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4
ATCGTCaaaGANNAATGANGgggTA[/mw_shl_code]
如果你计算机太烂了,对hg19基因组搞不定,就用这个测试数据也行。

对这个测试数据,我的处理如下:[mw_shl_code=perl,true]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
[/mw_shl_code]
如果是完整的perl代码:[mw_shl_code=perl,true]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"  

}[/mw_shl_code]
结果如下;
>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基因组下载方法有二十多种,我们统一用下面这个:
[mw_shl_code=shell,true]nohup wget http://hgdownload.cse.ucsc.edu/g ... Zips/chromFa.tar.gz &

tar zvfx chromFa.tar.gz

cat *.fa > hg19.fa

rm chr*.fa[/mw_shl_code]




上一篇:生信编程直播第五题:根据GTF画基因的多个转录本结构
下一篇:生信编程直播第三题:hg38每条染色体基因,转录本的分布
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

2

主题

34

帖子

777

积分

高级会员

Rank: 4

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

[mw_shl_code=applescript,true]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')[/mw_shl_code]

由于电脑问题,试了一下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

代码如下:
[mw_shl_code=python,true]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')[/mw_shl_code]
(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文件
脚本如下
[mw_shl_code=shell,true]tar zvfx chromFa.tar.gz   
cat *.fa > hg19.fa  
rm chr*.fa  
less hg19.fa [/mw_shl_code]按照老师的方法对python算法进行改良
改良后的代码如下:
[mw_shl_code=python,true]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))[/mw_shl_code]
保存代码为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可视化处理
[mw_shl_code=python,true]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+.35,s=str(round(cg_list)))
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()[/mw_shl_code]






回复 支持 1 反对 0

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

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

201-python-无名

自己的思考解答:学习了第一题之后,感觉python编程处理生物问题挺有意思的。第二题直接对长达3G的genome操作,有点吓人啊。幸好群主给了一个思路,先用小一点的文件测试下,输出结果为A、T、G、C、N的个数。我的思路:初始化一个字典,逐行读入文件并判断,如果是注释行,将染色体名加入字典并初始化一个空字符串,将注释行下的序列赋值给该染色体,文件读完后,对字典进行统计。用测试的小文件还是ok的,但是对基因组操作就一直卡在哪里。先把我的代码放在这里,学习老师的代码后再好好总结下。
我的代码如下:
[mw_shl_code=python,true]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])))
[/mw_shl_code]测试文件结果
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的视频:李老师讲解了两种方法。其中之一与我的方法类似,但是点拨了一下我们,即不能把所有的序列存入字典,这样占内存很大,可以逐个序列打印,我理解的逐个序列是按染色体来打印,因而结合老师的代码调整了一下,调整代码如下:

[mw_shl_code=python,true]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)
[/mw_shl_code]
测试文件结果:
>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的代码
[mw_shl_code=python,true]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)        
[/mw_shl_code]

运行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。
代码如下:
[mw_shl_code=python,true]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))[/mw_shl_code]
运行结果:
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

[mw_shl_code=python,true]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()
[/mw_shl_code]
回复 支持 3 反对 0

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-11 16:53:35 | 显示全部楼层
最近看书,正好看到Biostrings这个包,发现可以解决这个问题!就试了一下,挺好用!把使用过程贴出来供大家学习!(ps:R语言处于学习阶段,暂时还不会编程。)
[mw_shl_code=applescript,true]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[/mw_shl_code]

QQ截图20170111164954.png
QQ截图20170111165015.png
人生若只如初见!
回复 支持 3 反对 0

使用道具 举报

0

主题

8

帖子

179

积分

注册会员

Rank: 2

积分
179
发表于 2017-1-15 02:41:14 | 显示全部楼层
028-Python-Ryu
解题思路:
1.读取fasta文件,fasta文件每个不同的序列开头都有一个“>”号,因此可作为不同染色体的标记
2.统计序列长度使用len(),统计N、G、C个数使用.count()
3.一开始时考虑像李恒的readfq那样,先将序列存入一个变量,最后再统计,但是考虑需要内存以及速度较慢,最后直接每行处理,比较节省内存

完成代码后,再看老师的视频,最后再修改了一下自己的代码,主要是格式化输出以及计算GC含量时取出N的数目再计算。老师介绍的pysam很不错,又学习到新的知识!赞!

明天再做一下统计蛋白编码区的题目。

代码如下:
[mw_shl_code=python,true]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")))[/mw_shl_code]

运行时间在我的笔记本上大约2分钟左右
2017-01-15 02-22-02屏幕截图.png
最后结果如下:
[mw_shl_code=text,true]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
[/mw_shl_code]

代码和结果已上传Github,欢迎小伙伴们与我讨论,谢谢。
回复 支持 2 反对 0

使用道具 举报

1

主题

15

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-5 21:25:21 | 显示全部楼层
perl:043
[mw_shl_code=perl,true]#!/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";
}
[/mw_shl_code]

无标题.png
回复 支持 0 反对 1

使用道具 举报

0

主题

3

帖子

99

积分

注册会员

Rank: 2

积分
99
发表于 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确实很好用。
[mw_shl_code=python,true]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)))[/mw_shl_code]
回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

133

积分

注册会员

Rank: 2

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

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-1-12 23:03:45 | 显示全部楼层
本帖最后由 lilith098 于 2017-1-15 19:22 编辑

python+137
状态:还没有看视频,根据基础内容和第一题学到的内容尝试解答。

解题思路:
1 下载hg19的数据,因计算机配置原因,只选择chrY.fa 进行统计;
2 删除序列的空格
3 忽略开头字符“>”的行;
4 将所有列输出到一个字符串中;
5 针对字符串完成计数和输出。

[mw_shl_code=python,true]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)[/mw_shl_code]

输出报错:
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
分步输出:

print(seq)

print(seq)

计数正常显示

计数正常显示

seqLen = len(seq)报错

seqLen = len(seq)报错


结果:seqLen = len(seq) 此步骤报错,报错内容:TypeError: 'int' object is not callable

求教各位,这个地方报错的原因是什么?我实在没有看出来。

另外,本次作业只针对一个chrY.fa,没有考虑到题目中test.fa的类型。下一步将尝试处理此问题。

20170115
学习了超长的视频,代码如下:
[mw_shl_code=applescript,true]#!/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))
[/mw_shl_code]

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位浮点百分数的输出

问题:
为什么我的输出不是按照正常的顺序输出?




回复 支持 1 反对 0

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-1-15 14:52:19 | 显示全部楼层
本帖最后由 anlan 于 2017-1-15 15:01 编辑

034-perl+shell

test.fasta为测试数据

perl代码如下:
[mw_shl_code=perl,true]#!/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;[/mw_shl_code]
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

使用道具 举报

7

主题

32

帖子

216

积分

中级会员

Rank: 3Rank: 3

积分
216
发表于 2017-1-6 19:13:45 | 显示全部楼层
本帖最后由 babytong 于 2017-1-8 00:50 编辑

025 perl+python

[mw_shl_code=python,true]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)[/mw_shl_code]
测试文本的结果如下
ans2.png
欢迎关注https://www.jianshu.com/u/4f5e357a6212
回复 支持 反对

使用道具 举报

1

主题

41

帖子

287

积分

中级会员

Rank: 3Rank: 3

积分
287
发表于 2017-1-7 15:03:20 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-7 15:07 编辑

下载了1号染色体的数据进行一个测试,我的电脑h19估计跑不动,但没想到的是chr1一个文件也跑了很久。
解题思路:
1,提取fasta中的序列
2,计算GC和N的含量[mw_shl_code=python,true]#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'))
[/mw_shl_code]
回复 支持 反对

使用道具 举报

7

主题

32

帖子

216

积分

中级会员

Rank: 3Rank: 3

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

已更新代码与结果
欢迎关注https://www.jianshu.com/u/4f5e357a6212
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-1-8 22:01:49 | 显示全部楼层
139-python-佘璇
一开始也是想着楼上兄弟的做法,后来想着麻烦,看了下貌似挺简单的:
第二题 统计碱基个数.png
回复 支持 反对

使用道具 举报

11

主题

52

帖子

280

积分

中级会员

Rank: 3Rank: 3

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

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

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-9 21:05:45 | 显示全部楼层
本帖最后由 y461650833y 于 2017-1-9 21:07 编辑

最近一直在看perl小骆驼书,看这个题目的时候,顺便也上网找了一些相关资料,不过对哈希那里,有时候还是有点懵!不过还是模仿着写了出来,不过也是来来回回改了好几次!
[mw_shl_code=perl,true]#!/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";


}[/mw_shl_code]

QQ截图20170109210109.png
人生若只如初见!
回复 支持 反对

使用道具 举报

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函数对字典进行排序
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

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

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

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-4-4 17:12 , Processed in 0.057711 second(s), 35 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.