搜索
楼主: Jimmy

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

  [复制链接]

0

主题

13

帖子

93

积分

注册会员

Rank: 2

积分
93
发表于 2017-1-7 11:25:03 | 显示全部楼层
本帖最后由 tingting 于 2017-1-22 15:17 编辑

121-R/python-婷

1.22日mark一下,时隔两周,终于开始刷题交作业了。R还在补基础,争取下周交R的作业。今天先交python的作业。

解题思路:
本题目就是要将网上下载的CCDS.current.txt文件中倒数第二列中的起始位点和终止位点取出,加和所有(终止位点-起始位点)则可以计算出外显子的总长度。本题目要注意跳过第一行、每行最后有换行符需要去除、部分行没有外显子信息需要去除、部分外显子有重叠区域需要去除。

# coding: utf-8
# In[3]:
import os
import re
from operator import itemgetter
os.chdir('F:')

# In[ ]:
exonLength = 0
overlapExons = {}
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)

这一道题几乎抄袭老师的视频,下一道题争取先想想、试试再看视频。

最后还要感谢老师手把手教我入门,jupyter很好用,一步一步看视频装的、用的,谢谢!


回复 支持 反对

使用道具 举报

0

主题

6

帖子

53

积分

注册会员

Rank: 2

积分
53
发表于 2017-1-7 12:32:19 | 显示全部楼层
122-R-江
我也是完全0基础啊,群里共享的R的书还没买到呢,作业起码一两周去了,我尽力而为了。
回复 支持 反对

使用道具 举报

1

主题

16

帖子

152

积分

注册会员

Rank: 2

积分
152
发表于 2017-1-7 13:10:07 | 显示全部楼层
139 - Python - 佘璇

一开始我也想用字典除重,感觉有点麻烦,后来看到有个兄弟说道set,豁然开朗,用了两种正则匹配,发现findall效率并不比前面的低,难道是数据量太小的缘故?晚上有时间再用shell写下,感觉shell的话也很简单啊

PS:值得一提的是应该是length=end - start +1才对吧?写完我才发现的,所以就懒得再去截图了

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

1

主题

3

帖子

90

积分

注册会员

Rank: 2

积分
90
发表于 2017-1-7 13:54:16 | 显示全部楼层
我有一个问题,大家是不是没有考虑部分重叠的情况?

本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

4

主题

8

帖子

237

积分

中级会员

Rank: 3Rank: 3

积分
237
发表于 2017-1-7 14:27:26 | 显示全部楼层
本帖最后由 Davey 于 2017-1-7 18:09 编辑

081-perl+python    生信菜鸟刚接触perl和python编程语言有些基础用法还看不太懂,只能照着群主的代码敲了几遍运行测试一下,对这段代码有了基本的认识,但关于hash一块的用法还是有点看不懂,正努力查找小骆驼学习hash的用法。关于python的用法目前还不懂,只能先看下书学学基础,希望能早日摸些门道跟着各位大神好好学习
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
######格式化输出获取外显子的位置信息
open IN,"./CCDS.current.txt" or die $!;
open OUT,">getexonloc.txt" or die $!;
while(<IN>){
        chomp;
        next if /^#/;
        my @array=split;
        /\[(.*?)\]/;
        my @temp=split(/,/,$1);
        foreach(@temp){
                $_=~s/-/\t/;
                print "$array[2]\t$array[6]\t$array[0]\t$_\n";
                }
        }
close IN;
close OUT;

######对获取的外显子位子信息计算外显子的总长度
open loc,"getexonloc.txt" or die $!;
my (@F,%hash,$tmp);
while(<loc>){
        chomp;
        @F=split;
        foreach(($F[3]-$ARGV[0])..($F[4]+$ARGV[0])){
            $hash{"$F[2]:$_"}=1;
                }
        }
close loc;
$tmp++ foreach keys %hash;
print "$tmp\n";





回复 支持 反对

使用道具 举报

0

主题

5

帖子

115

积分

注册会员

Rank: 2

积分
115
发表于 2017-1-7 15:22:24 | 显示全部楼层
视频,看得很模糊哎。
这个R语言视频,apply等函数等看不懂。
惭愧,菜鸟。
回复 支持 反对

使用道具 举报

1

主题

4

帖子

160

积分

注册会员

Rank: 2

积分
160
发表于 2017-1-7 15:58:17 | 显示全部楼层
011-perl-早早
最近毕业论文实验好忙,估计这段时间跟不上了,但是我还是要为信仰充值。
回复 支持 反对

使用道具 举报

0

主题

3

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2017-1-7 17:08:40 来自手机 | 显示全部楼层
016-python+R
回复 支持 反对

使用道具 举报

0

主题

3

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2017-1-7 17:30:10 | 显示全部楼层
016-Python+R
目前还是小白,0基础,没怎么接触过编程,好多都不懂。Python敲了前两章,还有些不是很理解。看了视频觉得好厉害,自己下载了CCDS.txt文件,试了试群主以及其他人的代码,自己扣了很长时间,没运行出来。目前正在努力的理解去看代码。还得多练习多敲代码继续努力。
回复 支持 反对

使用道具 举报

0

主题

2

帖子

30

积分

新手上路

Rank: 1

积分
30
发表于 2017-1-7 17:36:36 | 显示全部楼层
本帖最后由 068-R 于 2017-1-8 21:24 编辑

学员编号:068-R

首先我们下载最新CCDs ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt
由于是0基础,R语言刚刚起步学,只能COPY 楼主的代码,其中几段语句并不是很懂,中年人表示压力好大~


[AppleScript] 纯文本查看 复制代码
rm(list=ls())
#按照老师的习惯清空数据

a=read.table("F:/coding/CCDS.current.txt",sep='\t',stringsAsFactors = F,header = T)  
#选择本地电脑目录的CCDs文件,以空格作为分间隔符,非factor模式,保留header

tmp <- apply(a[1:100,], 1, function(gene){# 取前100行数据分析调试 
 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, 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]) )# y[3]为坐标起点,y[4]为终止坐标,历编
        tmp_pos <<- c(tmp_pos,paste0(y[2],"-",i))
    })
    
  }
})
length(tmp_pos) # 计算exon的长度
length(unique(tmp_pos)) #计算去重后的exon的长度

由于电脑比较慢,运行许久还没跑出来
但中间出现
”In data.frame(gene = gene[3], chr = gene[1], start = start, end = end) :短变数里找到了行名,因此不用说是报错信息 不知道是什么原因

忘了说一点,打开CCDS.current.txt,删去第一行的#,不然读取数据表头会有问题
另外运行了x2yline的代码,速度好快!






回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-21 01:14 , Processed in 0.034412 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.