搜索
楼主: Jimmy

生信编程直播第一题:人类基因组的外显子区域到底有多...

  [复制链接]

0

主题

4

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-1-6 10:02:59 | 显示全部楼层

非常谢谢你的回复啊
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-6 10:09:13 | 显示全部楼层
huigezi056 发表于 2017-1-6 10:02
非常谢谢你的回复啊

这个在群主的视频里面 他说了
人生若只如初见!
回复 支持 反对

使用道具 举报

7

主题

32

帖子

216

积分

中级会员

Rank: 3Rank: 3

积分
216
发表于 2017-1-6 11:01:14 | 显示全部楼层
learnyoung 发表于 2017-1-5 22:02
059 python+R-冷洋
试了一些方法,一开始用存入字典的方法来去重,发现效率有点低,我的笔记本CORE_M的处理 ...

你的python代码中漏了f.close()
推荐使用with 来打开文件,因为它在用完文件后会自动释放掉所占的内存
欢迎关注https://www.jianshu.com/u/4f5e357a6212
回复 支持 反对

使用道具 举报

7

主题

32

帖子

216

积分

中级会员

Rank: 3Rank: 3

积分
216
发表于 2017-1-6 12:05:17 | 显示全部楼层
本帖最后由 babytong 于 2017-5-2 16:00 编辑

025 perl +python

[Python] 纯文本查看 复制代码
import re
sum=0
list_all=[]

with open ('D:/sxjns/question1/CCDS.current.txt') as f:
    for line in f:
        if line.startswith('#'):continue
        lists=line.rstrip().split('\t')
        exons=re.findall('[0-9]+-[0-9]+',lists[-2])
        [list_all.append(exon) for exon in exons if not exon in list_all]  
    for i in list_all:
        exon_one=i.split('-')
        sum+=int(exon_one[1])-int(exon_one[0])
    print(sum)


1.with打开文件优于传统方法打开
2.事先判断存在性,再append 就不需要进行去重
3.以上代码运算结果为36027374 以上代码执行相当慢
所以做出如下修改
[Python] 纯文本查看 复制代码
import os
import re
os.chdir("D:")
sum=0
list_all=[]
with open('CCDS.current.txt','r') as f:
    for line in f :
        if line.startswith('#'):continue
        lists=line.rstrip().split('\t')
        exons=re.findall('[0-9]+-[0-9]+',lists[-2])
        for exon in exons:
            list_all.append(exon) 
        
    list_all=set(list_all)
    for i in list_all:
        exon_one=i.split('-')
        sum+=int(exon_one[1])-int(exon_one[0])
    print(sum)

025 perl +python


欢迎关注https://www.jianshu.com/u/4f5e357a6212
回复 支持 2 反对 0

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-1-6 14:29:59 | 显示全部楼层
本帖最后由 chapman 于 2017-1-6 14:56 编辑

050-Python-晁帅
[Python] 纯文本查看 复制代码
import refrom collections import OrderedDict

sum = 0
overlap = OrderedDict()

fi = open('CCDS.current.txt','rt')
for each_line in fi:
    if each_line.startswith('#'):
        continue
    each_line = each_line.rstrip()#去掉句尾换行符
    line = each_line.split('\t')#按换行符分开
    if line[-2] == '-':
        continue
    line_ex = re.sub("\[|\]"," ",line[-2])#去掉[]


    target_lines = line_ex.split(",")#去掉,
    for target_line in target_lines:
        start = int(target_line.split('-')[0])
        end = int(target_line.split('-')[1])
        coordinate = line[0] + ':'+ target_line
        if coordinate not in overlap.keys():
            overlap[coordinate] = 1    
            sum += end - start
fi.close()
print(sum)
    



1.不怎么会用with打开,个人感觉with打开有点牛逼;2.如果py文件和目标文件不在一个目录下,还是要用os.chdir()的;
3.最终结果是:36054296,和大家的不一样,不知道为啥?

回复 支持 反对

使用道具 举报

0

主题

3

帖子

28

积分

新手上路

Rank: 1

积分
28
发表于 2017-1-6 16:04:14 | 显示全部楼层
本帖最后由 end2end 于 2017-1-6 16:20 编辑

I used some tricks like eval and sum,but I'm not sure if 925941-926012 contains 926012 or not?




[Python] 纯文本查看 复制代码
import pandas
import re
table=pandas.read_table("CCDS.current.txt")
string="".join(table.values[:,-2])
l=re.findall("\d+-\d+",string)
num=-sum([eval(i) for i in set(l)])

image.png

by 075-python-哈哈哈
回复 支持 1 反对 0

使用道具 举报

0

主题

4

帖子

49

积分

新手上路

Rank: 1

积分
49
发表于 2017-1-6 17:44:36 | 显示全部楼层
049-R+python+shell-刘王宓
还算顺利,得出了视频中同样的结果

import re
import os
from collections import OrderedDict
from operator import itemgetter

os.chdir('D:/')

exonLength = 0
overlapExons = OrderedDict()
with open('CCDS.current.txt','rt') as f:
    for line in f:
        if line.startswith('#'):
            continue
        line = line.rstrip()
        lst = line.split('\t')
        if lst[-2]  == '-':
            continue
        lst[-2] = re.sub('\[|\]','',lst[-2])
        exons = lst[-2].split(', ')
        for exon in exons:
            start = int(exon.split('-')[0])
            end = int(exon.split('-')[1])
            coordinate = lst[0]+':'+exon
            if coordinate not in overlapExons.keys():
                overlapExons[coordinate] = 1
                exonLength += end-start

print (exonLength)

36048075
回复 支持 反对

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-6 20:24:52 | 显示全部楼层
learnyoung 发表于 2017-1-5 22:44
把外显子存入list,然后元组set一下去重,10s之内结果就出来了

谢谢!

修改了一下代码:

结果是一样的,但是2秒就出来结果了:)

修改了两个地方:
1、将strip()替换成了re.sub()
2、引用了set()去重

[Python] 纯文本查看 复制代码
import re

path = '/Users/sunshine/practice/practice1-计算人类基因组coding区的总长度/'

f_raw = open(path+'CCDS.current.txt','r')
f_raw_lines = f_raw.readlines()[1:]

cds_length = 0
position = []

for line in f_raw_lines:
    line = line.strip('\n').split('\t')
    cds_locations = line[9]
    if cds_locations != ('-'):
        cds_locations = re.sub('\[|\]','',cds_locations).split(', ')
        for cds_location in cds_locations:
            position.append(cds_location)
            
position = list(set(position))
for cds_location in position:
    cds_location = cds_location.split('-')
    cds_length += (int(cds_location[1])-int(cds_location[0])+1)

print (cds_length)
f_raw.close()
                               
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-1-6 22:57:38 | 显示全部楼层
learnyoung 发表于 2017-1-5 22:02
059 python+R-冷洋
试了一些方法,一开始用存入字典的方法来去重,发现效率有点低,我的笔记本CORE_M的处理 ...

厉害了 我的哥
回复 支持 反对

使用道具 举报

2

主题

34

帖子

777

积分

高级会员

Rank: 4

积分
777
发表于 2017-1-7 08:59:39 | 显示全部楼层
本帖最后由 x2yline 于 2017-2-18 10:29 编辑

077 R-python- shell  x2yline
1.我的解题思路:




2.python实现
(1)jupyter编辑器太强大了,非常好用,但是没有查看当前变量的功能,所以最终还是选择spyder作为python编写平台(有shift+enter键相当于Rstudiod的ctr+r键,也有查看当前已有变量数值的功能
(2)关于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来表示换行。

(3)代码实现
[Python] 纯文本查看 复制代码
import pandas as pd
import numpy as np
file = 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

3.R语言实现

第一次是未考虑外显子重叠情况写的R语言代码如下:
运行计算时间为:14.74084 secs
最后运行结果为:36048075

Rplot02.png
代码如下:
[AppleScript] 纯文本查看 复制代码
setwd('E:\\r\\biotrainee_demo1')
t1 <- Sys.time()
directory = 'CCDS.current.txt'

data <- read.table(directory, sep='\t',
                   stringsAsFactors=F, header=T)[c(1,10)]

get_gene <-function(data_item){
  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){
  exon <- unique(strsplit(gene,", ")[[1]])
}

get_length <- function(exon){
  loc <- strsplit(exon,"-")[[1]]
  a <- as.numeric(loc[2])-as.numeric(loc[1])
  if (a==0){
    print(loc)
  }
  a
}


exon_length = 0
exon_length_items = NULL
for (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)
  print('************************')
  print(paste('染色体', 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
}

# 耗时长度
difftime(Sys.time(), t1, units = 'secs')

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

经过差不多一下午的思考和调试(大多数时间都在改写代码,因为很多种方法R语言都跑不动

最后终于也做出来了,速度也还可以:
运行计算时间为:18.71507 secs
最后运行结果为: 34889110
运行时间:2017-01-07 18:06:38 CST
系统配置:
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


代码如下:
[AppleScript] 纯文本查看 复制代码
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 = 0
for (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 <- exon_length + dim(base_loc)[1] * dim(base_loc)[2]
}

# 耗时长度
difftime(Sys.time(), t1, units = 'secs')

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












回复 支持 2 反对 0

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-20 21:08 , Processed in 0.040980 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.