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

题目 人类基因组外显子区域长度

学员:x2yline

具体题目详情请参考生信技能树论坛

题目数据来源为:ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt

解题思路(比较适合R语言)如下

R语言实现

第一次完成的代码是未考虑外显子overlap情况(只去重了完全相同的外显子)写的

运行计算时间:14.74084 secs

最后运行结果:36048075

第一版代码如下:

setwd('E:\\r\\biotrainee_demo\\class1')#修改工作路径t1 <- Sys.time()#把程序运行之前的时间赋值给t1directory = 'CCDS.current.txt'#把文件名赋值给directorydata <- read.table(directory, sep='\t',                   stringsAsFactors=F, header=T)[c(1,10)]#读取数据并提取出第一和第十列get_gene <-function(data_item){  # 该函数用于apply执行  # 输入的数据为仅含原始数据第1列和第10列的dataframe  # 用apply函数执行后输出的数据为每个基因外显子的坐标,  # 一个基因的所有外显子以逗号分隔组成一个string,所有基因的string组成一个vector  # 用apply函数执行后,最后格式为c('111-112, 115-135, 125-138', '254-258',...)  if (!data_item[2] =='-'){    exon_ranges <- data_item[2]    exon_ranges <- substr(exon_ranges, start=2, stop=nchar(exon_ranges)-1)# 去除首尾的中括号符号  }}get_exon <- function(gene){  # 输入的数据为c('111-112, 115-135', '125-138', '254-258,...')  # 把i号染色体上的所有外显子后在一起,并去除完全相同的外显子  # 输出的数据为c('111-112','115-135', '125-138', '254-258', ...)  exon <- unique(strsplit(gene,", ")[[1]])}get_length <- function(exon){  # 输入的数据为lapply(c('111-112','115-135', '125-138', '254-258', ...),fun)  # 输出结果为左右两坐标之差+1即外显子的长度  loc <- strsplit(exon,"-")[[1]]  a <- as.numeric(loc[2])-as.numeric(loc[1]) +1 #每个外显子的碱基数目  a}exon_length = 0exon_length_items = NULLfor (i in unique(data[,1])){  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')  exon_i <-  get_exon(gene_i)  exon_i_length <- sapply(exon_i, get_length)  exon_length <- exon_length + sum(exon_i_length)  exon_length_items <- c(exon_i_length, exon_length_items)  names(exon_length_items)[1:length(exon_i_length)] <- i}hist(exon_length_items,xlim=c(0,500),breaks = 20000,      main='Distribution of exon length', xlab='exon length')difftime(Sys.time(), t1, units = 'secs')# 计算执行完成后时间与t1的间隔print(paste('all exons length is',exon_length))

第二版代码如下

setwd('E:\\r\\biotrainee_demo1')t1 <- Sys.time()directory = 'CCDS.current.txt'# 读取数据并提取第1列和第10列data <- read.table(directory, sep='\t',                   stringsAsFactors=F, header=T)[c(1,10)]get_gene <-function(data_item){  # 用apply执行该函数  # 输入的数据为仅含原始数据第1列和第10列的dataframe  # 输出的数据为c('111-112, 115-135, 125-138', '254-258',...)  if (!data_item[2] =='-'){    exon_ranges <- data_item[2]    exon_ranges <- substr(exon_ranges, start=2, stop=nchar(exon_ranges)-1)  }}get_exon <- function(gene){  # 输入的数据为c('111-112, 115-135, 125-138, 254-258,...')  # 输出的数据为c('111-112','115-135', '125-138', '254-258', ...)  exon <- unique(strsplit(gene,", ")[[1]])# 注:strsplit的输出结果为列表}get_length <- function(exon){  # 输入的数据为lapply(c('111-112','115-135', '125-138', '254-258', ...),fun)  # 输出结果为两坐标值和左右两坐标之差  loc <- strsplit(exon,"-")[[1]]  a <- c(as.numeric(loc[1]), as.numeric(loc[2])-as.numeric(loc[1]), as.numeric(loc[2]))  #if (a==0){  #print(loc)  #}  a}exon_length = NULLfor (i in unique(data[,1])){  # paste 函数把i号染色体的所有外显子的坐标合并为一个character对象  # gene_i的格式为'111-112, 115-135, 125-138, 254-258,...'  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')  # exon_i的格式为c('111-112','115-135', '125-138', '254-258', ...)  exon_i <-  lapply(get_exon(gene_i), get_length)  mat <- matrix(unlist(exon_i), ncol=3, byrow = T)  #mat <- mat[order(mat[,2], decreasing = F),]  #mat <- mat[order(mat[,1], decreasing = F),]  # 使用matrix 是因为vector太长会报错  #R memory management / cannot allocate vector of size n MB  base_loc <- matrix(unique(unlist(apply(mat, 1, function(x) c(x[1]:x[3])))))  exon_length <- c(exon_length , dim(base_loc)[1] * dim(base_loc)[2])}# 耗时长度difftime(Sys.time(), t1, units = 'secs')chrs <- unique(data[,1])barplot(exon_length,names.arg=chrs,xlab='Chromosomes',ylab='length of exons')print(paste('all exons length is',sum(exon_length)))

python实现

  • jupyter编辑器太强大了,非常好用,但是没有查看当前变量的功能,所以最终还是选择spyder作为python编写平台(有shift+enter键相当于Rstudiod的ctr+r键,也有查看当前已有变量数值的功能)
  • 关于open(file, 'rt')的解释

w,r,wt,rt都是python里面文件操作的模式。w是写模式,r是读模式。

t是windows平台特有的所谓text mode(文本模式),区别在于会自动识别windows平台的换行符。

类Unix平台的换行符是\n,而windows平台用的是\r\n两个ASCII字符来表示换行,python内部采用的是\n来表示换行符。

rt模式下,python在读取文本时会自动把\r\n转换成\n. wt模式下,Python写文件时会用\r\n来表示换行。

python代码实现(第一次写的与老师的代码大致相同,用for循环即可,不推荐用以下的方法做)

import pandas as pdimport numpy as npfile = r'E:\r\biotrainee_demo1\CCDS.current.txt'def calculate_exon(file):  data = pd.read_csv(file, sep='\t',\    usecols=[0,9])#data.loc[1:10,:]#  data[0:3]#  data.iloc[1:3]#  data.iloc[3]  all_length = 0  for i in data.iloc[:,0].unique():    # get the data of chrosome i    # iloc[row_vector,col_vect]    # iloc[row_vector]    data_i = data.loc[data.iloc[:,0] == i]    type(data_i)    type(data_i.iloc[:,1])    # remove the '[]' in column2    data_j = data_i.iloc[:,1].apply(lambda x: x[1:-1])    data_p = data_j.apply(lambda x: x.split(', '))    data_g = data_p.apply(lambda x: pd.Series(x))    # 把nan填充为 0-0    data_f = np.array(data_g.fillna('0-0'))    # 去除重复的外显子    data_f = np.unique(data_f.reshape((data_f.shape[0]*data_f.shape[1], 1)))    data_f = pd.DataFrame(data_f)    data_m = data_f.apply(lambda x: \      x.apply(lambda y: (y.split('-')[0])))    data_n = data_f.apply(lambda x: \      x.apply(lambda y: (y.split('-')[-1])))    # pd.to_numeric can only apply to a 1-d array    data_mi = data_m.apply(lambda x: pd.to_numeric(x, downcast='float'))    data_ni = data_n.apply(lambda x: pd.to_numeric(x, downcast='float'))    all_length += (data_ni - data_mi).sum().sum()  return(all_length)length = calculate_exon(file)print(length) 

运算速度有点慢,因为是临时学的pandas和numpy,很多步骤还没有优化

未去重overlap结果为:36046283

编程感悟

由于开始R是没有基础的,用通过R包swirl学习了一下lapply,apply和sapply函数的使用,对于迭代数目比较多的循环来说,R语言的for循环效率远远不如apply系列函数,应该尽量避免for循环处理,而python的for循环运算速度较快,可以使用for循环处理一下比较大的数据。

本文编辑:思考问题的熊

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-03-15

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏小樱的经验随笔

深入理解树状数组

树状数组(Binary Indexed Tree(BIT), Fenwick Tree)是一个查询和修改复杂度都为log(n)的数据结构。主要用于查询任意两位之...

2527
来自专栏包子铺里聊IT

[Google最新面试题] continental divider

给一个矩阵,其中0代表海洋,其他数字代表高度, 秉着水往低处流的原则,求出能够 流向任意海洋的点。 比如说 0 0 0 1 2 3 0 0 1 2 2 4 3 ...

2554
来自专栏Java架构沉思录

如何优雅地过滤敏感词

敏感词过滤功能在很多地方都会用到,理论上在Web应用中,只要涉及用户输入的地方,都需要进行文本校验,如:XSS校验、SQL注入检验、敏感词过滤等。今天着重讲讲如...

411
来自专栏xingoo, 一个梦想做发明家的程序员

汇编语言 手记4

简单的汇编指令 ? CPU执行后,寄存器中的数据改变为如下: ? CPU访问内存单元时要给出内存单元的地址。所有的内存单元构成的存储空间是一个一维的线性空间。 ...

1715
来自专栏数据分析

char varchar nchar nvarcharar到底有多大区别

首先说明下,ASP.NET MVC系列还在龟速翻译中。 工作好多年,基础知识甚是薄弱,决定以后在coding(cv操作)的时候尽量多google下,然后总结下来...

2656
来自专栏挖掘大数据

处理海量数据的10种常见方法

本文将介绍10种处理海量数据问题的常见方法,也可以说是对海量数据的处理方法进行一个简单的总结,希望对你有帮助。

20210
来自专栏SeanCheney的专栏

《图解算法》总结第1章 算法简介第2章 选择排序第3章 递归第4章 快速排序第5章 散列表第6章 广度优先搜索第7章 狄克斯特拉算法第8章 贪婪算法第9章 动态规划

Grokking Algorithms: An illustrated guide for programmers and other curious peop...

3559
来自专栏蜉蝣禅修之道

LeetCode之Jump Game II

1324
来自专栏大数据挖掘DT机器学习

【知识】SAS数据分析完整笔记(3)

SAS学习笔记(3):SAS一般高级语言 本篇SAS读书笔记主要介绍SAS一般高级语言,主要内容包括赋值语句、输出语句、分支机构、循环结构、数组以及函数等六个...

3489
来自专栏算法channel

常用排序算法代码兑现

主要推送关于对算法的思考以及应用的消息。坚信学会如何思考一个算法比单纯地掌握100个知识点重要100倍。本着严谨和准确的态度,目标是撰写实用和启发性的文章,欢迎...

33811

扫描关注云+社区