生信编程直播课程优秀学员作业展示2
时间:2022-05-03
本文章向大家介绍生信编程直播课程优秀学员作业展示2,主要内容包括学员:x2yline、R语言实现(太卡,高能报警)、python3第一种实现方法(运行速度较快,但没有3快)、基本概念、基础应用、原理机制和需要注意的事项等,并结合实例形式分析了其使用技巧,希望通过本文能帮助到大家理解应用这部分内容。
题目:hg19基因组序列的一些探究
学员:x2yline
具体题目详情请参考生信技能树论坛
数据来源:http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz 下载.gz数据后解压
R语言实现(太卡,高能报警)
代码地址:https://raw.githubusercontent.com/x2yline/courseranotes/master/myscript/class2/fastafileGCandN.R
代码内容:
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')
由于电脑问题,试了一下1号染色体,电脑卡住了,于是又试了一下Y染色体,跑出来结果如下:
耗时:41.44945 secs
电脑配置信息:
- 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
python3第一种实现方法(运行速度较快,但没有3快)
数据来源: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
数据下载时间:2017-01-10 23:08
运行消耗时间:309 seconds
未优化速度的代码如下
import osimport timebegin = time.time()os.chdir(r'F:tmpchromFa')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) + ' secondsn')
shell +python3(最快)
先使用shell脚本把所有chromFa.tar.gz 中的所有.fa文件合并为一个hg19.fa文件
脚本如下:
tar zvfx chromFa.tar.gzcat *.fa > hg19.farm chr*.faless hg19.fa
按照老师的方法对python算法进行改良
改良后的代码如下:
代码地址:
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 = 'chromsometGC_contenttN_contenttLengthtNtAtTtCtGn' 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 ('nnSorry! 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:tmphg19.fa
部分输出结果如下
使用python进一步进行可视化处理
代码如下:
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 = 'chromsometGC_contenttN_contenttLengthtNtAtTtCtGn' 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:genomechromFahg19.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()
本文编辑:思考问题的熊
- 一条SQL语句的执行计划变化探究(r10笔记第3天)
- tensorflow(一)windows 10 python3.6安装tensorflow1.4与基本概念解读
- 基于AgileEAS.NET SOA 中间件领域模型数据器快速打造自己的代码生成器
- Java基础-day07-代码题-自定义数据类型;ArrayList集合
- 一条报警信息的快速处理和分析(r9笔记第99天)
- 【Go 语言社区】解析Go语言编程中的struct结构
- centos+scala2.11.4+hadoop2.3+spark1.3.1环境搭建
- 【Go 语言社区】Golang 语言获取本机逻辑CPU数量的方法
- Data Guard搭建困境突围(一)(r10笔记第17天)
- Java基础-day07-知识点相关题-自定义数据类型;ArrayList
- windows10 tensorflow(二)原理实战之回归分析,深度学习框架(梯度下降法求解回归参数)
- 本人为巨杉数据库(开源NoSQL)写的C#驱动,支持Linq,全部开源,已提交github
- 最近的几个技术问题总结和答疑(九)(r10笔记第16天)
- AgileEAS.NET SOA中间件平台更新日志 2015-04-28
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 浅析动态切换数据源的原理(接上篇)
- SpringBoot源码解析(十二)- Autowired是如何注入的
- 项目要实现多数据源动态切换,咋搞?
- 这一次,带你全面了解锁机制!
- GitHub标星1w+超牛的微服务项目,开发脚手架
- Redis中hash、set、zset的底层数据结构原理
- Redis中string、list的底层数据结构原理
- Redis中字符串的表示
- Redis分布式锁背后的原理
- 解析Transformer模型
- 这5个常问的Redis面试题你答得出来吗?(详细剖析)
- 性能最佳实践:MongoDB索引
- Python基本数据类型-list-tuple-dict-set
- 深度学习应用的服务端部署
- MongoDB中的CURD操作