生信编程直播课程优秀学员作业展示2

setwd('E:\\r\\biotrainee_demo\\class 2')# 读入数据t1 <- Sys.time()df <- read.csv('chr1.fa', 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所在的位置+1index_df1 <- data.frame(rbind(matrix(index_df[-1,1]),dim(df)[1]+1))# 把index_start和index_end存入data.frameindex_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 <- strsplit(y, split=character(0))[[1]]  N <- length(y[y =='N'])/length(y)  GC <- length(y[y =='G' | y == 'C'])/(length(y)-length(y[y =='N']))  c(N,GC)})# 把行名改为N和GC并转秩rownames(result) = c('N','GC')result <- t(result)# 取结果前几行head(result)difftime(Sys.time(), t1, units = 'secs')

import osimport timebegin = 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 += nall_info = 'chr,GC_ratio,N_ratio,Length,N,A,T,C,G'for i in range(len(chromsome)):    data = '\n'+str(chromsome[i]) +',' + "%.5f"%((content[i][-1]+content[i][-2])/sum(content[i][1:])) +','  + "%.5f" %(content[i][0]/(sum(content[i]))) +','  +str((sum(content[i]))) +','  +str((content[i][0])) + ','  +str(content[i][1])+',' +str(content[i][2])+','  +str(content[i][3])+','  +str(content[i][4])     all_info += datawith open('hg19_analysis.csv','w') as f:    f.write(all_info)print('Time using:'+ str(time.time() - begin) + ' seconds\n')

tar zvfx chromFa.tar.gzcat *.fa > hg19.farm chr*.faless hg19.fa

import osimport timeimport reimport sysfrom collections import OrderedDictstart = 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()        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'):].upper()                        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_analysisdef 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_contentif sys.argv:    result = OrderedDict()    for f in sys.argv:        done = 0        f= f.strip(''''"''')        if f.count('.') != 1 or 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))

python fasta_atcgn_summary.py   F:\tmp\hg19.fa

import osimport timeimport reimport sysfrom collections import OrderedDictstart = 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_analysisdef 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_contentfile_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 pltplt.bar(left = range(25), height = cg_list, color='k')for i in range(len(cg_list)):    plt.text( x=i-0.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:\hg19_gc.png',dpi=600)plt.show()