搜索
楼主: Jimmy

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

  [复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-1-7 17:56:09 | 显示全部楼层
002-Python-R-shell

python 自己写还吃不消,但是调试了东哥,以及其他朋友的几种方法。因为代码不是自己原创的就不粘贴了。

收获:
1.仔细研究了正则表达式
2.re.sub()函数
3. python 文件的读取,以及路径(path)问题
4. file() 函数

还差的太远,需要多多努力~!
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-1-7 18:03:07 | 显示全部楼层
023-R 不用谈了
R学习刚刚开始,处理这个问题真的很有难度,不过根据视频代码能跑出来。这几天在恶补R的基础知识,R in action 看了三章,收获很大,希望尽快理解视频中代码的意思。
首先,我也是下载了文件名是CCDS.current.txt文件,大概9.1M。
a<-read.table("CCDS.current.txt",sep='\t',stringsAsFactors = F,head=T)
head(a)后,我发现一个问题,第一行也就是列标签没显示!下载不同版本的CCDS都是一样,我觉得很奇怪,勉强打开后发现我下载的CCDS首行是以#开头,这可能导致R不识别,删除后成功显示!
然后以下是运行结果:
tmp<-apply(a[1:100,],1,function(gene){
    #gene=a[3,]
    x=gene[10]#column10外显子所在序列
    if(grepl('\\]',x)){#判断X中是否有这一符号,如果有就用正则替换掉
      x=sub('\\[','',x)
    x=sub('\\]','',x)
    tmp<-strsplit(as.character(x),',')[[1]]#我们先从逗号开始分割
    start<-as.numeric(unlist(lapply(tmp,function(y){#取开始点
      strsplit(as.character(y),'-')[[1]][1]})))
    end<-as.numeric(unlist(lapply(tmp,function(y){#取结束点
      strsplit(as.character(y),'-')[[1]][2]})))
    gene_d<-data.frame(#构建数据框
  gene=gene[3],chr=gene[1],start=start,end=end)
   return(gene_d)#返回数据框
    }})
------
  tmp_pos=c()#构造一个空向量
lapply(tmp[1:10],function(x){
  #print(x)
  if(!is.null(x)){
    apply(x,1,function(y){
      #print(y)
      for ( i in as.numeric(y[3]):as.numeric(y[4]) )
    tmp_pos<<-c(tmp_pos,paste0(y[2],"-",i))
    }
      )
  }
})#运行后结果返回

$`1`
NULL

$`2`
NULL

$`3`
NULL

$`4`
NULL

$`5`
NULL

$`6`
NULL

$`7`
NULL

$`8`
NULL

$`9`
NULL

$`10`
NULL
length(unique(tmp_pos))
[1] 10087。#结果显示是这样,我需要理解以下
回复 支持 反对

使用道具 举报

0

主题

6

帖子

125

积分

注册会员

Rank: 2

积分
125
发表于 2017-1-7 18:05:41 | 显示全部楼层
osho 发表于 2017-1-7 18:03
023-R 不用谈了
R学习刚刚开始,处理这个问题真的很有难度,不过根据视频代码能跑出来。这几天在恶补R的基 ...

不过是不是不太对,少这么多??
回复 支持 反对

使用道具 举报

2

主题

34

帖子

773

积分

高级会员

Rank: 4

积分
773
发表于 2017-1-7 18:13:43 | 显示全部楼层
068-R 发表于 2017-1-7 17:36
学员编号:068-R

首先我们下载最新CCDs ftp://ftp.ncbi.nlm.nih.gov/pub/C ... n/CCDS.20160908.txt

我的代码又改了一下
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-7 19:15:17 | 显示全部楼层
osho 发表于 2017-1-7 18:05
不过是不是不太对,少这么多??

这个是10个基因外显子的长度 直播的老师 没有算全部的 太慢
人生若只如初见!
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-7 19:22:12 | 显示全部楼层
本帖最后由 y461650833y 于 2017-1-8 16:30 编辑

027-R+Perl-游戏
今天晚上看了直播 怎么用R完成第一题 顺手就把代码都打了下来。最近在看R语言实战,才把前面7章看完,因此对于直播老师用的很多函数,还不是很熟悉,不过这个代码,还是可以看懂的!
setwd("E:/perl")
getwd()
[1] "E:/perl"
rm(list = ls())
dir()
[1] "CCDS.current.txt"
a=read.table("CCDS.current.txt",sep = '\t',stringsAsFactors = F,header = T)##记得把下载文件第一行开始的#号去掉 要不会读不出来的##

tmp <- apply(a, 1, function(gene){
  #gene=a[3,]
  x=gene[10]
  if(grepl('\\]',x)){
    x=sub('\\[','',x)
    x=sub('\\]','',x)
    tmp <- strsplit(as.character(x),',')[[1]]
    start <- as.numeric(unlist(lapply(tmp, function(y){
      strsplit(as.character(y),'-')[[1]][1]
    })))
    end <- as.numeric(unlist(lapply(tmp, function(y){
      strsplit(as.character(y),'-')[[1]][2]
    })))
    gene_d <- data.frame(gene=gene[3],
                         chr=gene[1],
                         start=start,
                         end=end
    )

    return(gene_d)
  }
})
##44M##
tmp_pos=c()
lapply(tmp, function(x){
  if (!is.null(x)){
    apply(x, 1, function(y){
      for (i in as.numeric(y[3]) :as.numeric(y[4]))
        tmp_pos <<- c(tmp_pos,paste0(y[2],"-",i))
    })
  }
})
length(tmp_pos)
length(unique(tmp_pos))
##结果暂时没运行出来,太慢了。。因为算的是全部的外显子长度!等运算出来了 再贴图!!ps:可能是笔记本太差,最后死机了,好尴尬,没看到结果!

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
人生若只如初见!
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-7 19:48:38 | 显示全部楼层
x2yline 发表于 2017-1-7 18:13
我的代码又改了一下

你的代码确实快!怎么做到的。。。
人生若只如初见!
回复 支持 反对

使用道具 举报

0

主题

3

帖子

57

积分

注册会员

Rank: 2

积分
57
发表于 2017-1-7 19:59:01 | 显示全部楼层
009-卉卉自己几乎没有基础,还在学习基础语法中。。。最近期末考试,感谢解答。

wget  ftp://ftp.ncbi.nlm.nih.gov/pub/C ... n/CCDS.20160908.txt  
vim tem.pl  
输入代码如下:
[Perl] 纯文本查看 复制代码
while(<>){
          next if /^#/;                                  #遇到#开头跳过
          @F=split;
          /\[(.*?)\]/;                                      #匹配 [ ]中内容
          @tmp=split/,/,$1;                          #这里为什么是$1?
          foreach(@tmp){
                  $_=~s/-/\t/;                          # -替换成空格
                  print "$F[2]\t$F[6]\t$F[0]\t$_\n";
          }
}

保存退出
perl tem.pl CCDS.20160908.txt


wc tmp.txt
外显子35万个

perl -alne'{$tmp+=($F[4]-$F[3])}END{print $tmp}' tmp.txt
长度 57M

perl -alne'{foreach ($F[3]..$F[4]) {$h{$_}=1} }END{$tmp++ foreach keys %h;print $tmp}'  tmp.txt


很遗憾了。。。。等两天换个服务器再试试






本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

1

主题

6

帖子

434

积分

中级会员

Rank: 3Rank: 3

积分
434
发表于 2017-1-7 20:20:59 | 显示全部楼层
041-shell+python
目前还没有什么编程基础,不过正在看Python的书和廖雪峰Python网站教程学习,跟习题一视频学习了一下基本理解了老师的思路和每行代码的作用,目前正处于照葫芦画瓢阶段,也利用老师的代码实际操作成功了,记录一下第一次课程的学习收获吧。
1.文件的读取;
2.条件语句和循环语句;
3.提取目标列后,还要对外显子分割,终止位置坐标-起始位置坐标后求和,就是怎么去重复的那块代码还没有特别明白,我再理解理解;
4.小目标:坚持自学编程基础,特别是时间充裕时更不能放松,学习生信编程能有群里的人指引还是非常幸运的。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-1-7 20:45:25 | 显示全部楼层
perl+python-137

解题思路:
1.     提取文件中每一行的cds location;
2.     去除格式,将location中起始坐标和终止坐标分别存为整数;
3.     行:排除:第一行:表头;其他行:内容为“-”;
4.     去除并列关系的外显子重复;
5.     用终止坐标减去起始坐标,结果累加得到外显子长度。

代码记录:
# coding: utf-8

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

os.chdir('D:/Programing/study/biotrainee/test 1')

exonLen = 0
overlapExons = OrderedDict()
with open('CCDS.current.txt','rt') as f:
    for line in f:  #every line
        if line.startswith('#'):
            continue
        line = line.rstrip() # removed the LF in the end of the line
        list = line.split('\t')  # spilt the line by tab
        if list[-2] == '-':
            continue
        list[-2] = re.sub('\[|\]', ' ', list[-2]) # Remove the [] in locations
        exons = list[-2].split(', ') # spilt the exon list by ",",exons is a list, like [1082892-1082986,1084352-1084382]
        for exon in exons:  #exon:1082892-1082986
            lenth = exon.split('-')    #lenth is a list,[1082892,1082986]
            start = int(lenth[0])     #1082892
            end = int(lenth [1])     #1082986
            cooridinate = list[0]+':'+exon  #record the chromosome and exons locations
            if cooridinate not in overlapExons.keys():
                overlapExons[cooridinate] = 1
                exonLen += end-start
print (exonLen)

输出结果:
36061455
小记:
只看完了python简明教程,所以对于视频中各种模块的用法不太清楚,需要尽快学习目前用到的模块。re:regular expression,正则表达式模块。OrderDict():有序的字典(数据类型)
问题:
1 operator 中itemgetter的目的。另外,我在pycharm+python2.7中,貌似没有这个模块,再加入去重复那一段代码后,程序一直运行不出来。暂时不知道是不是这个问题。
2 f在这里是作为一个table的对象吗?
3 会不会出现1082892-1082986,1082893-1082998或1082983-1082900等,部分重叠或完全包含的情况?


回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 15:42 , Processed in 0.033511 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.