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

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

学员：x2yline

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

题目数据来源为：

解题思路（比较适合R语言）如下

R语言实现

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

运行计算时间：14.74084 secs

最后运行结果：36048075

第一版代码如下：

1. setwd('E:\\r\\biotrainee_demo\\class1')#修改工作路径

2. t1 <- Sys.time()#把程序运行之前的时间赋值给t1

3. directory = 'CCDS.current.txt'#把文件名赋值给directory

4. data <- read.table(directory, sep='\t',

5.                   stringsAsFactors=F, header=T)[c(1,10)]#读取数据并提取出第一和第十列

6. get_gene <-function(data_item){

7.  # 该函数用于apply执行

8.  # 输入的数据为仅含原始数据第1列和第10列的dataframe

9.  # 用apply函数执行后输出的数据为每个基因外显子的坐标，

10.  # 一个基因的所有外显子以逗号分隔组成一个string，所有基因的string组成一个vector

11.  # 用apply函数执行后，最后格式为c('111-112, 115-135, 125-138', '254-258',...)

12.  if (!data_item[2] =='-'){

13.    exon_ranges <- data_item[2]

14.    exon_ranges <- substr(exon_ranges, start=2, stop=nchar(exon_ranges)-1)# 去除首尾的中括号符号

15.  }

16. }

17. get_exon <- function(gene){

18.  # 输入的数据为c('111-112, 115-135', '125-138', '254-258,...')

19.  # 把i号染色体上的所有外显子后在一起，并去除完全相同的外显子

20.  # 输出的数据为c('111-112','115-135', '125-138', '254-258', ...)

21.  exon <- unique(strsplit(gene,", ")[[1]])

22. }

23. get_length <- function(exon){

24.  # 输入的数据为lapply(c('111-112','115-135', '125-138', '254-258', ...),fun)

25.  # 输出结果为左右两坐标之差+1即外显子的长度

26.  loc <- strsplit(exon,"-")[[1]]

27.  a <- as.numeric(loc[2])-as.numeric(loc[1]) +1 #每个外显子的碱基数目

28.  a

29. }

30. exon_length = 0

31. exon_length_items = NULL

32. for (i in unique(data[,1])){

33.  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')

34.  exon_i <-  get_exon(gene_i)

35.  exon_i_length <- sapply(exon_i, get_length)

36.  exon_length <- exon_length + sum(exon_i_length)

37.  exon_length_items <- c(exon_i_length, exon_length_items)

38.  names(exon_length_items)[1:length(exon_i_length)] <- i

39. }

40. hist(exon_length_items,xlim=c(0,500),breaks = 20000,

41.     main='Distribution of exon length', xlab='exon length')

42. difftime(Sys.time(), t1, units = 'secs')# 计算执行完成后时间与t1的间隔

43. print(paste('all exons length is',exon_length))

第二版代码如下

1. setwd('E:\\r\\biotrainee_demo1')

2. t1 <- Sys.time()

3. directory = 'CCDS.current.txt'

4. # 读取数据并提取第1列和第10列

5. data <- read.table(directory, sep='\t',

6.                   stringsAsFactors=F, header=T)[c(1,10)]

7. get_gene <-function(data_item){

8.  # 用apply执行该函数

9.  # 输入的数据为仅含原始数据第1列和第10列的dataframe

10.  # 输出的数据为c('111-112, 115-135, 125-138', '254-258',...)

11.  if (!data_item[2] =='-'){

12.    exon_ranges <- data_item[2]

13.    exon_ranges <- substr(exon_ranges, start=2, stop=nchar(exon_ranges)-1)

14.  }

15. }

16. get_exon <- function(gene){

17.  # 输入的数据为c('111-112, 115-135, 125-138, 254-258,...')

18.  # 输出的数据为c('111-112','115-135', '125-138', '254-258', ...)

19.  exon <- unique(strsplit(gene,", ")[[1]])# 注：strsplit的输出结果为列表

20. }

21. get_length <- function(exon){

22.  # 输入的数据为lapply(c('111-112','115-135', '125-138', '254-258', ...),fun)

23.  # 输出结果为两坐标值和左右两坐标之差

24.  loc <- strsplit(exon,"-")[[1]]

25.  a <- c(as.numeric(loc[1]), as.numeric(loc[2])-as.numeric(loc[1]), as.numeric(loc[2]))

26.  #if (a==0){

27.  #print(loc)

28.  #}

29.  a

30. }

31. exon_length = NULL

32. for (i in unique(data[,1])){

33.  # paste 函数把i号染色体的所有外显子的坐标合并为一个character对象

34.  # gene_i的格式为'111-112, 115-135, 125-138, 254-258,...'

35.  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')

36.  # exon_i的格式为c('111-112','115-135', '125-138', '254-258', ...)

37.  exon_i <-  lapply(get_exon(gene_i), get_length)

38.  mat <- matrix(unlist(exon_i), ncol=3, byrow = T)

39.  #mat <- mat[order(mat[,2], decreasing = F),]

40.  #mat <- mat[order(mat[,1], decreasing = F),]

41.  # 使用matrix 是因为vector太长会报错

42.  #R memory management / cannot allocate vector of size n MB

43.  base_loc <- matrix(unique(unlist(apply(mat, 1, function(x) c(x[1]:x[3])))))

44.  exon_length <- c(exon_length , dim(base_loc)[1] * dim(base_loc)[2])

45. }

46. # 耗时长度

47. difftime(Sys.time(), t1, units = 'secs')

48. chrs <- unique(data[,1])

49. barplot(exon_length,names.arg=chrs,xlab='Chromosomes',ylab='length of exons')

50. 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循环即可，不推荐用以下的方法做）

1. import pandas as pd

2. import numpy as np

3. file = r'E:\r\biotrainee_demo1\CCDS.current.txt'

4. def calculate_exon(file):

5.  data = pd.read_csv(file, sep='\t',\

6.    usecols=[0,9])

7. #data.loc[1:10,:]

8. #  data[0:3]

9. #  data.iloc[1:3]

10. #  data.iloc[3]

11.  all_length = 0

12.  for i in data.iloc[:,0].unique():

13.    # get the data of chrosome i

14.    # iloc[row_vector,col_vect]

15.    # iloc[row_vector]

16.    data_i = data.loc[data.iloc[:,0] == i]

17.    type(data_i)

18.    type(data_i.iloc[:,1])

19.    # remove the '[]' in column2

20.    data_j = data_i.iloc[:,1].apply(lambda x: x[1:-1])

21.    data_p = data_j.apply(lambda x: x.split(', '))

22.    data_g = data_p.apply(lambda x: pd.Series(x))

23.    # 把nan填充为 0-0

24.    data_f = np.array(data_g.fillna('0-0'))

25.    # 去除重复的外显子

26.    data_f = np.unique(data_f.reshape((data_f.shape[0]*data_f.shape[1], 1)))

27.    data_f = pd.DataFrame(data_f)

28.    data_m = data_f.apply(lambda x: \

29.      x.apply(lambda y: (y.split('-')[0])))

30.    data_n = data_f.apply(lambda x: \

31.      x.apply(lambda y: (y.split('-')[-1])))

32.    # pd.to_numeric can only apply to a 1-d array

33.    data_mi = data_m.apply(lambda x: pd.to_numeric(x, downcast='float'))

34.    data_ni = data_n.apply(lambda x: pd.to_numeric(x, downcast='float'))

35.    all_length += (data_ni - data_mi).sum().sum()

36.  return(all_length)

37. length = calculate_exon(file)

38. print(length)

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

未去重overlap结果为：36046283

编程感悟

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

---

本文编辑：思考问题的熊