搜索
楼主: Jimmy

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

  [复制链接]

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-1-5 13:39:50 | 显示全部楼层
[Python] 纯文本查看 复制代码
import os
import re

os.chdir("F:/生信菜鸟团作业")
f1=open("CCDS.current.txt")

exonlenghth=0
line=f1.readlines()

for i in line:
    if i.startswith("#"):
        continue
    i=i.rstrip()
    filt=re.split(r'\t',i)
    if filt[-2]=="-":
        continue
    z=re.sub("\[|\]"," ",filt[-2])
    exons=z.split(",")
    for exon in exons:
        start=int(exon.split("-")[0])
        end=int(exon.split("-")[1])
        exonlenghth+=end-start
print (exonlenghth)
回复 支持 反对

使用道具 举报

0

主题

4

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-1-5 21:07:07 | 显示全部楼层
next if /^#/;这句话到底是什么意思呢?036问
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-5 21:09:25 | 显示全部楼层
huigezi056 发表于 2017-1-5 21:07
next if /^#/;这句话到底是什么意思呢?036问

跳过#开头的行
人生若只如初见!
回复 支持 反对

使用道具 举报

1

主题

41

帖子

285

积分

中级会员

Rank: 3Rank: 3

积分
285
发表于 2017-1-5 22:02:17 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-8 00:49 编辑

059 python+R-冷洋
试了一些方法,一开始用存入字典的方法来去重,发现效率有点低,我的笔记本CORE_M的处理器跑了好久才出来结果(大概是去细胞房处理了一批细胞的时间)。所以想能不能提高效率一点,每一次通过判断exon在不在字典里,不在还要写进去来去重,太慢了(对于 我这个老爷本来说是这样,等的好焦急)。在python中set就是一个没有重复的集合,我把用正则 [0-9]+-[0-9]+匹配到外显子存入一个list中,然后将list 经过set化就去重了,验证了一下速度非常快,比写入字典的方法要快很多解题思路:
1,正则[0-9]+-[0-9]+直接逐行匹配CDS_location
2,将匹配到的CDS_location存入一个list中
3,set(list)去重复
4,首尾相减,算长度
方法一:
[Python] 纯文本查看 复制代码
#coding=utf-8
import re
f=open('CCDS.current.txt','r')
sum=0
list=[]
dit={}
for line in f:
    if line.startswith('#'):
        continue
    line=line.rstrip()
    exons=re.findall('[0-9]+-[0-9]+',line)#匹配像00000-00000这种格式的数据,检查了一下可以将所有的外显子匹配到
    for exon in exons:
        list.append(exon)#将每一个外显子存入list中
st=set(list)#将列表set去重
for i  in st:
    exon_one=i.split('-')#还可以用find函数找到-的位置然后用切片操作找到起始点和终止点
    exon_start=exon_one[0]
    exon_end = exon_one[1]
    sum+=int(exon_end)-int(exon_start)
print sum
。结果为36027064

方法二:
[Python] 纯文本查看 复制代码
#coding=utf-8
import re
f=open('CCDS.current.txt','r')
g=open('NCCDS.txt','w')
sum=0
dit={}
for line in f:
    if line.startswith('#'):
        continue
    line=line.rstrip()
    exons=re.findall('[0-9]+-[0-9]+',line)
    for exon in exons:
        exon_one=exon.split('-')
        exon_start=exon_one[0]
        exon_end=exon_one[1]
        if exon not in dit.keys():
            dit[exon]=1
        sum+=int(exon_end)-int(exon_start)
print sum

[Python] 纯文本查看 复制代码
#coding=utf-8
import re
f=open('CCDS.current.txt','rt')
sum=0
lt=[]
dic={}
for line in f:
    line=line.rstrip()#去掉每一行最后的换行符
    if line.startswith('#'):
        continue
    lst=line.split('\t')#以空格将
    if lst[-2]=='-':
        continue
    exons=re.sub("\[|\]"," ",lst[-2])#去除中括号
    exons=exons.split(',')
    for exon in exons:
        lt.append(exon)
st=set(lt)
for i in st:
    exon_one=i.split('-')
    exon_start = exon_one[0]
    exon_end = exon_one[1]
    sum += int(exon_end) - int(exon_start)
print sum







回复 支持 5 反对 0

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-1-5 22:16:51 | 显示全部楼层
python-60

基于CCDS.current.tx文件,去除了重复的coding区,计算得到coding区总长度为36237316,结果较视频教材中的结果(36048075)多,我想其原因应该是我在每次计算的时候都加了1(如,第一个coding区,视频中为71(926012-925941),可我则为72。

可是计算花的时间有点多,花了快11分钟,不知道是哪里耽误了计算,求指导。

代码如下:

[Python] 纯文本查看 复制代码
path = '/Users/sunshine/practice/'

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

cds_length = 0
position = []

for line in f_lines:
    line = line.strip('\n').split('\t')
    cds_locations = line[9]
    if cds_locations != ('-'):
        cds_locations = cds_locations.strip("'").strip('[').strip(']').split(', ')
        for cds_location in cds_locations:
            if cds_location not in position:
                position.append(cds_location)
                cds_location = cds_location.split('-')
                cds_length += (int(cds_location[1])-int(cds_location[0])+1)

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

使用道具 举报

1

主题

41

帖子

285

积分

中级会员

Rank: 3Rank: 3

积分
285
发表于 2017-1-5 22:44:52 | 显示全部楼层
Teanna2017 发表于 2017-1-5 22:16
python-60

基于CCDS.current.tx文件,去除了重复的coding区,计算得到coding区总长度为36237316,结果较 ...

把外显子存入list,然后元组set一下去重,10s之内结果就出来了
回复 支持 反对

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

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

dit[exon]=1兄弟,这点为什么要等于1啊
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-1-5 23:28:30 | 显示全部楼层
034-perl-shell

文件是:CCDS.20160908.txt

代码如下:

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $in = "CCDS.20160908.txt";
my $out = "CCDS.20160908_exon.txt";

open my $fh, $in or die "Cannot open $in";
open my $fh_out, ">$out" or die "Cannot write $out";
readline $fh;
while (<$fh>){
	chomp;
	my @array = split /\t/, $_;
	if ($_ =~ /\[(.*?)\]/){
		my @tmp = split /,/, $1;
		map{
			$_ =~ s/-/\t/;
			print $fh_out "$array[2]\t$array[6]\t$array[0]\t$_\n";
		}@tmp;
	}
}
close $fh;
close $fh_out;

my $in2 = "CCDS.20160908_exon.txt";
my %hash;
open $fh, $in2 or die "Cannot open $in2";
while (<$fh>){
	chomp;
	my @array = split /\t/, $_;
	foreach ($array[3]..$array[4]){
		$hash{"$array[0]:$_"} = 1;
	}
}
close $fh;

my $length = keys %hash;
print "$length\n";


结果为:35011112
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-1-5 23:31:43 | 显示全部楼层
本帖最后由 y461650833y 于 2017-2-14 10:55 编辑

这次的作业,只能看看大家的作业了
刚学perl不久,看的小骆驼,刚把前面三章看完!
在看群主以及其他童鞋的代码的时候,好多函数啥的还不是很明白,为了看懂,把小骆驼书是一顿翻啊
特别是哈希求外显子长度之和哪里,是真的不懂
运行了几个大家的作业代码,发现可以运行的通。
自己还需要努力学习027-R+perl-游戏 perl小白一个,求解救!##趁着今天有空,把视频里面的代码敲了一遍##
#!/usr/bin/perl -w
while(<>){
          next if /^#/;                                 
          @F=split;
          /\[(.*?)\]/;                                      
          @tmp=split/,/,$1;                          
          foreach(@tmp){
                  $_=~s/-/\t/;                        
                  print "$F[2]\t$F[6]\t$F[0]\t$_\n";
          }
}
## perl zy1_1.pl CCDS.current.txt > zuobiao.txt
## perl -alne'{$tmp+=($F[4]-$F[3])}END{print $tmp}' zuobiao.txt
57746136
##perl -alne'{foreach ($F[3]..$F[4]) {$h{$_}=1} }END{$tmp++ foreach keys %h;print $tmp}'  zuobiao.txt


##由于用的是自己的笔记本,提示内存溢出##
#!/usr/bin/perl -w
while(<>){
  chomp;
  @zuobiao=split;
  foreach ($zuobiao[3]..$zuobiao[4]){
    $hash{"$zuobiao[2] : $_"}=1;
  }
}
$sum++ foreach keys %hash;
print "$sum\n";

#perl zy1_2.pl zuobiao.txt #内存溢出 也是醉了###

##对于利用哈希求长度,还不是很熟悉,需要继续学习,其他方面,倒是都理解了##


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

使用道具 举报

1

主题

41

帖子

285

积分

中级会员

Rank: 3Rank: 3

积分
285
发表于 2017-1-6 09:24:54 | 显示全部楼层
017中大普外 发表于 2017-1-5 23:10
dit[exon]=1兄弟,这点为什么要等于1啊

去重的时候我们把已经计算过长度的外显子存入到字典dic中,下次一个外显子要计算长度时先判断他在不在字典中个,如果在就不计算长度(去重),不在就计算。创建字典的时候字典里每一个元素必须是以key:value的形式存在,这里dic[exon]=1就是为了满足这个条件,1不重要,起一个建立字典的辅助作用,写成dic[exon]=2也没关系
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.