搜索
查看: 3481|回复: 178

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

  [复制链接]

360

主题

650

帖子

2275

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
2275
发表于 2016-8-27 21:20:39 | 显示全部楼层 |阅读模式
首先,我们要知道外显子的坐标文件,参考:http://www.bio-info-trainee.com/1101.html
可以看到根据NCBI的CCDS参考,统计到的外显子全才是34729283bp,也就是约35M,根很多文献里面说的都一样!同理,我统计了一下外显子的侧翼长度
54692160 外显子加上前后50bp
73066288  外显子加上前后100bp
90362533  外显子加上前后150bp

而hg19版本基因组里面有着entrez gene ID号的基因是23056个基因,所以我接下来探究一下这些基因的信息!
我们首先看看基因与基因之间的交叉情况
其中有12454266bp的位点,是多个外显子共有的,可能是一个基因的多个外显子,或者是不同基因的外显子
我们主要看看不同基因的交叉情况
不同基因交叉情况共516种,共涉及498种基因!
ZNF341   NFS1
ZNF397   ZSCAN30
我们来看看,如果最不擅长的R来做,会是怎么样的呢?
[Python] 纯文本查看 复制代码
rm(list=ls())

a=read.table(choose.files(),sep = '\t',stringsAsFactors = F,header = T)# 选择你下的CCDs文件

tmp <- apply(a[1:100,], 1, function(gene){# 取前100行数据分析调试
  # gene=a[3,]
  x=gene[10]# Column10 外显子坐标位置列
  if(grepl('\\]',x)){#判断x中是否存在有]这样的符号,如果有就利用正则替换掉。
    x=sub('\\[','',x)
    x=sub('\\]','',x)
    # 这个时候得到的对象还是像这样的“880073-880179, 880436-880525……”
    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){ #取前10个list文件计算调试
 # 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的长度


用perl我只需要两个单行命令,三五分钟就搞定了,可是在R里面,写这个代码花了我足足半个小时,因为我的R学的很差,而且,最后压根没有做出来!!
但是这个代码很有借鉴意义,你们好好听视频吧!
里面所有的R的apply系列代码,请务必拆开来看!!!

语言,本质上都是一样的,用R语言也可以编程制作一个QQ聊天软件。
实现一个需求,在R里面各种函数都可以,也可以用各种封装好的包,这就是我对R语言的实验,用包!







上一篇:人类hg19这个参考基因组里面各个染色体长度是多少?
下一篇:人类参考基因组中所有基因对应的转录本分布情况如何呢?
回复

使用道具 举报

0

主题

6

帖子

75

积分

注册会员

Rank: 2

积分
75
发表于 2017-1-8 00:03:06 | 显示全部楼层
028-Python-Ryu

1. 什么是外显子?什么是cds?
搜索后在biostars看到
CDS: "A contiguous sequence which begins with, and includes, a start codon and ends with, and includes, a stop codon."

http://www.sequenceontology.org/ ... cvs/term/SO:0000316

Exon: "A region of the transcript sequence within a gene which is not removed from the primary RNA transcript by RNA splicing."

http://www.sequenceontology.org/ ... cvs/term/SO:0000147
那么个人观点是,外显子是形成mRNA后剪接剩下的部分,包括UTR区与CDS,而CDS则是真正编码蛋白的Coding Sequence;
2.此次作业虽然是统计人类基因组的外显子区域,但是还是统计CDS部分,不统计外显子;
3.去重不去除在染色体上overlap的部分,只去除完全一样的部分;

1.首先,根据老师的做法是,从NCBI上下载CCDS的最新版本文件,ftp://ftp.ncbi.nlm.nih.gov/pub/C ... an/CCDS.current.txt
代码如下
[Python] 纯文本查看 复制代码
import csv

ncbi_file="CCDS.20160908.txt"
with open(ncbi_file,'r')as f1:
    file=csv.reader(f1,delimiter="\t")
    next(file)
    sum=0
    exon_dict={} # dict is faster than list
    for record in file:
        if record[9] != "-":
            chr=record[0]
            exon_list=record[9].lstrip("[").rstrip("]").split(", ")
            # print(exon_list)
            for range in exon_list:
                exon=chr+":"+range
                # print(exon)
                if exon not in exon_dict:
                    exon_dict[exon]=""
                    exon_start=int(range.split("-")[0])
                    exon_end=int(range.split("-")[1])
                    # print(exon_start)
                    sum+=exon_end-exon_start
print(sum)
运行结果
36M,与老师的相符,虽然结果一致,但是如果仅仅是重复老师做的,那就没有意思了;其实在看老师的视频前,这道题我想到的方法是去UCSC或者Ensembl下载所有基因的区域再进行统计,老师从NCBI下载算是教会我NCBI还有这样的公共库,下面我就从自己的方法去验证结果:

2.从UCSC下载所有CDS的region来验证,下载方法是使用ucsc的table browser,链接在这里
下载后,文件打开,是这样的,统计“Exon_Starts”和“Exon_Ends”即为每个基因的cds起始和结束位置
代码如下
[Python] 纯文本查看 复制代码
import csv

ucsc_file="ucsc_cds.txt"
with open(ucsc_file,'r')as f1:
    exon_dict={}
    sum=0
    file=csv.reader(f1,delimiter="\t")
    next(file)
    for record in file:
        chr=record[1]
        exon_start_list=record[5].rstrip(",").split(",")
        exon_end_list=record[6].rstrip(",").split(",")
        pos=0
        while pos < len(exon_start_list):
            exon=chr+":"+exon_start_list[pos]+"-"+exon_end_list[pos]
            if exon not in exon_dict:
                exon_dict[exon]=""
                sum+=(int(exon_end_list[pos])-int(exon_start_list[pos]))
            pos+=1
print(sum)
运行结果
33M,基本相符

3.从Ensembl下载所有的CDS进行验证,使用的是Ensembl的Biomart,链接在这里
Biomart需要点选的数据库为
从上图可以看出,如果统计“Exon Chr Start”“Exon Chr End”会把一些非编码的区域也统计进去,需要统计的是“Genomic Coding Start”和“Genomic Coding End”,当然,取巧的统计“CDS Length”这一列应该也可以,但是要去重,我还没有尝试,我试着比较统计Exon和CDS的两个结果
代码如下
[Python] 纯文本查看 复制代码
import csv

ensembl_file="mart_export.txt"
with open(ensembl_file,'r')as f1:
    file=csv.reader(f1,delimiter="\t")
    next(file)

    cds_sum=0
    exon_sum=0

    cds_dict={} # dict is faster than list
    exon_dict={}

    for record in file:
        chr=record[3]
        exon_start=record[5]
        exon_end=record[6]
        exon=chr+":"+exon_start+"-"+exon_end
        if exon not in exon_dict:
            exon_dict[exon]=""
            exon_sum+=(int(exon_end)-int(exon_start))

        if record[7] != "":
            cds_start=record[7]
            cds_end=record[8]
            cds=chr+":"+cds_start+"-"+cds_end
            if cds not in cds_dict:
                cds_dict[cds]=""
                cds_sum+=(int(cds_end)-int(cds_start))

print(cds_sum)
print(exon_sum)
运行结果,先输出CDS总长度再输出Exon总长度
可以看到,CDS结果是39M,基本符合预期。但是外显子统计结果则有117M,相差的有点大,估计是数据有些冗余,有基因的不同剪接和非编码的区域

欢迎各位同学指点,谢谢!


本帖子中包含更多资源

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

x
回复 支持 12 反对 0

使用道具 举报

0

主题

1

帖子

32

积分

新手上路

Rank: 1

积分
32
发表于 2017-1-17 21:31:12 | 显示全部楼层
105-R+Python
刚刚开始接触python ,此题以学习并理解老师的代码为主:
1.生物基本概念:
1)基因DNA分为编码区和非编码区,编码区包含内含子和外显子,一般非编码区具有基因表达的调控功能,如启动子在非编码区。编码区则转录为mRNA并最终翻译成蛋白质。
2)内含子和外显子都被转录到mRNA的前提hnRNA中,当hnRNA进行剪接变为成熟的mRNA时,内含子被切除,外显子被保留。实际上真正编码蛋白质的是外显子,而内含子则无编码功能。
3)CDS Sequencecoding for aminoacids in protein 蛋白质编码区 外显子并不一定编码蛋白质
2.Python
以研究学习别人的代码为主
1)Python模块
在Python中,一个.py文件就被称之为一个模块。
使用模块的最大好处提高了代码的可维护性,以及代码不必从零开始
模块的引用:模块名.函数名
2)代码中的3个模块
A.OS模块
OS 操作系统   与操作系统相关的功能 ,且不受平台限制
常用功能:
        • 显示当前使用平台 os.name()
        • 显示当前python 脚本工作路径 os.getcwd()
        • 删除一个文件  os.remove('filename')  
        • 重命名文件 os.rename("oldname","newname") 
B.CSV模块
csv模块是python的内置模块,直接import csv 即可
csv.reader()---读取csv文件数据
csv.writer()---写入csv文件数据
示例:1.improt csv  #读取csv
                 reader=csv.reader(open('test.csv','rb'))
                 for item in reader:
                        print line
               2.import csv  #写入csv
                  writer=csv.writer(open('test.csv','wb'))---此处必须加上参数'b',不然会出现各行现象
C.RE模块
Regular Expression 正则表达式
提供各种正则表达式的匹配操作,在文本解析、复杂字符串分析和信息提取时是一个非常有用的工具   
3)next()函数
 next()方法当一个文件被用作迭代器,典型例子是在一个循环中被使用,next()方法被反复调用。此方法返回下一个输入行,或引发StopIteration异常EOF时被命中。
for index in range(5):
  line = fo.next()
  print "Line No %d - %s" % (index, line)
4)with open(…) as myfile
对open的调用,会返回一个简单的文件对象,复制给myfile。可以使用一般的文件工具来使用myfile。此对象也支持with语句所使用的环境管理协议。在这个with语句执行后,环境管理机制保证由myfile所引起的文件对象会自动关闭,即使处理该文件时,for循环引发了异常也是如此。
                             


回复 支持 2 反对 0

使用道具 举报

0

主题

14

帖子

143

积分

注册会员

Rank: 2

积分
143
发表于 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


代码如下:
[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))












本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

4

主题

27

帖子

145

积分

注册会员

Rank: 2

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

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

025 perl +python
回复 支持 2 反对 0

使用道具 举报

1

主题

3

帖子

57

积分

注册会员

Rank: 2

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

本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

0

主题

1

帖子

18

积分

新手上路

Rank: 1

积分
18
发表于 2017-2-8 11:07:44 | 显示全部楼层
006-perl
第一课就看了好几遍,理解了老师的思路,但是在os下的操作尚未掌握,未能成功运行。在学习共享的Linux 和perl语言的基础知识资料。寒假已经补了一些Linux基本命令。希望自己能把基础学习好后完成作业
回复 支持 1 反对 0

使用道具 举报

360

主题

650

帖子

2275

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
2275
 楼主| 发表于 2017-1-17 09:06:22 | 显示全部楼层
朝曦曦 发表于 2017-1-17 05:31
大家是怎么上传代码的?

http://www.biotrainee.com/forum-62-1.html 请自己修改一下
回复 支持 1 反对 0

使用道具 举报

0

主题

1

帖子

16

积分

新手上路

Rank: 1

积分
16
发表于 2017-1-10 09:59:46 | 显示全部楼层
本帖最后由 tony2015116 于 2017-1-16 19:34 编辑

这几周忙着复习考试,所以没做作业,抱歉了视频看完了,表示压力山大,很多地方不明白
回复 支持 1 反对 0

使用道具 举报

0

主题

3

帖子

63

积分

版主

Rank: 7Rank: 7Rank: 7

积分
63
发表于 2017-1-9 10:24:17 | 显示全部楼层

拟南芥外显子长度

本帖最后由 fei0810 于 2017-1-9 10:35 编辑

本身是做植物表观的,所以这里没有再分析人的外显子区域,而是算了模式植物 Arabidopsis
计算使用TAIR 10 的gff 文件,并没有使用最新的Araport 11数据,下载地址 :https://www.arabidopsis.org/down ... IR10_GFF3_genes.gff
由于是GFF文件所以格式比较整齐,不需要做太多的预处理,就是一个行筛选、行加减、去重复最后列求和的问题,可以分别用Python 和 linux 自带的命令来进行一下尝试

[Python] 纯文本查看 复制代码
import os
import re
from operator import itemgetter
from collections import OrderedDict

os.chdir('I:python/homework1/')

exonLength = 0
overlapExons = OrderedDict()
with open('./TAIR10_GFF3_genes.gff','rt')as file:
    for line in file:
        lst = line.split("\t")
        if lst[2] == 'CDS':
            start = int(lst[3])
            end = int(lst[4])
            coordinate = lst[0]+':'+lst[3]+lst[4]
            if coordinate not in overlapExons.keys():
                    overlapExons[coordinate] = 1
                    exonLength += end-start

print (exonLength)


结果是34,686,354

如果使用Linux 的话看起来更简洁一些,写法也可以说各种各样。

[Bash shell] 纯文本查看 复制代码
cat TAIR10_GFF3_genes.gff |awk 'BEGIN{OFS="\t"}{if($3=="CDS") print $1,$3,$4,$5,$5-$4}'| sort -k1,1 -k3,3n|uniq |awk '{sum += $5};END {print sum}' 


结果是34,686,353

2000年拟南芥全基因组的Nature文献里给的数据是33,249,250

084-Python+R-熊

回复 支持 1 反对 0

使用道具 举报

1

主题

40

帖子

202

积分

中级会员

Rank: 3Rank: 3

积分
202
发表于 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







回复 支持 1 反对 0

使用道具 举报

0

主题

12

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2017-1-3 20:41:41 | 显示全部楼层
用perl写了一个,对于所注释的行做的不够好。

回复 支持 反对

使用道具 举报

0

主题

12

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2017-1-3 20:42:40 | 显示全部楼层
本帖最后由 btrainee 于 2017-1-3 20:43 编辑

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

open  IN ,"CCDS.current.txt" or die $!;
my $loc;
my %hash;
while (<IN>){
>---chomp;
>---next if /^#/;
>---$_=~/\[(.*?)\]/;
>---$loc.=$1.",";
}
close IN;
my @location=split /,/,$loc;
foreach my $tt (@location){
#>--next if $tt=="-";
>---my @tmp=split /-/,$tt;
>--->---foreach my $i ($tmp[0]..$tmp[1]){
>--->--->---$hash{$i}=1;
>--->---}
>---}
my $hash_size=keys %hash;
print "$hash_size\n";
回复 支持 反对

使用道具 举报

3

主题

19

帖子

180

积分

版主

Rank: 7Rank: 7Rank: 7

积分
180
发表于 2017-1-4 10:57:15 | 显示全部楼层
回复

使用道具 举报

1

主题

7

帖子

78

积分

注册会员

Rank: 2

积分
78
发表于 2017-1-4 16:01:10 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w

use strict;
my @temp;
my @nextdata;
my @nextnextdata;
my $zuobiao;
my $site;
my %hash;
my $sum=1;
open my $data,"CCDS.current.txt" or die $!;

while(<$data>){
chomp;
next if /^#/;
 @temp=split/\t/,$_;
if ($temp[9]=~/\[(.*)\]/){
 @nextdata=split/,\s/,$1;
foreach $zuobiao (@nextdata){
 @nextnextdata=split/-/,$zuobiao;
foreach  $site ($nextnextdata[0]..$nextnextdata[1])
{$hash{$site}=1}
}
}
}
foreach (keys %hash){$sum++;}
print $sum;
close $data;


计算结果31691416
perl-024
回复 支持 反对

使用道具 举报

0

主题

12

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2017-1-4 20:32:03 | 显示全部楼层
本帖最后由 btrainee 于 2017-1-4 20:35 编辑

048-perl+R+python
一开始写的版本有点问题,把所有重复的坐标全去除了,没有考虑到不同的染色体。重新写了个版本,使用最新下载的20170101CCDS.current.txt,如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;
open  IN ,"CCDS.current.txt" or die $!;
my $loc;
my %hash;
while (<IN>){
>---chomp;
>---next if /^#/;
>---my @F=split;
>---(my $match)=$_=~/\[(.*?)\]/;
>---$match=~s/(\d+)\-(\d+)/$F[0]\-$1\-$2/g;
>---$loc.=$match.",";
}
close IN;
my @location=split /,/,$loc;
foreach my $tt (@location){
#>--next if $tt=="-";
>---my @tmp=split /-/,$tt;
>--->---foreach my $i ($tmp[1]..$tmp[2]){
>--->--->---$hash{"$tmp[0]:$i"}=1;-
>--->---}
>---}
my $hash_size=keys %hash;
print "$hash_size\n";

运行结果为:35288998
然后又发现一个问题,把版主的perl,copy过来运行,发现结果不一致
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;
#############格式化输出
my (@F,$match,@L,%hash,$tmp);
open IN,"CCDS.current.txt" or die $!;
open OUT,">","CCDS.20170101.exon.txt" or die $!;
while (<IN>){
>---next if /^#/;
>---@F=split;
>---($match)=$_=~/\[(.*?)\]/;
>---my @tmp=split/,/,$match;
>---foreach (@tmp){
>--->---$_=~s/-/\t/g;
>--->---print OUT "$F[2]\t$F[6]\t$F[0]\t$_\n";
>---}
}
close IN;
close OUT;
#############算外显子总长度,或者加上侧翼50/100/150后的总长度
open FH,"CCDS.20170101.exon.txt" or die $!;
while (<FH>){
>---chomp;
>---@L=split;
>---foreach (($L[3]-$ARGV[0])..($L[4]+$ARGV[0])){
>--->---$hash{"$L[2]:$_"}=1;
>---}
}
close FH;
$tmp++ foreach keys %hash;
print "$tmp\n";

$ARGV[0]参数给为0,运行结果为:34889110
由于python还不太会,所以直接copy了视频的代码,自己运行了一下,神奇了,结果也不同,代码如下:
[Python] 纯文本查看 复制代码
import re
import os
from collections import OrderedDict
from operator import itemgetter

os.chdir('~/bob/btraine/1_exon_length')
exonlength = 0
overlapExons = OrderedDict()
with open ('CCDS.current.txt','rt') as f:
>---for line in f:
>--->---if line.startswith('#'):
>--->--->---continue
>--->---line = line.rstrip() #去掉右边的换行符空格TAB等,lstrip strip
>--->---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
一时有点晕,用的同一个文件,怎么计算结果都不一样,有点费解。
回复 支持 反对

使用道具 举报

3

主题

7

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2017-1-4 20:52:27 | 显示全部楼层
本帖最后由 猪兔子2号 于 2017-1-4 21:25 编辑

摩羯座:007

首先下载文件,linux 命令行下输入 wget -bc ftp://ftp.ncbi.nlm.nih.gov/pub/C ... n/CCDS.20160908.txt

然后写脚本,前面跟群主的一样

[Perl] 纯文本查看 复制代码
#命令行输入perl test.pl CCDS.20160908.txt >test.txt
while(<>){
               chomp;
               next if /^#/;
               @F=split("\t",$_);
                /\[(.*?)\]/;
                @tmp=split(",",$1);
          foreach (@tmp){ 
         @region=split("-",$_);
          print "$F[2]\t$F[6]\t$F[0]\t$region[0]\t$region[1]\n";                                   } 
              }#命令行输入perl test1.pl test.txt >test2.txt
$i=0
while(<>){
             @tmp=split("\t",$_);
             foreach ($tmp[3]..$tmp[4]) {
                  if($hash{"$tmp[2]:$_"}){print $hash{"$tmp[2]:$_"}."\t"$tmp[0]}
                  else{
                            $hash{"$tmp[2]:$_"}=$tmp[2]; $i++;}                                                                                       
         }                         
   }print $i;


这里不一样的地方时,命名一个变量,起始为零,判断该位点是否存在,不存在加一,输出所有外显子位点之和,并且输出那些基因具有重复位点两端的基因有哪些,因为有些基因不知道为什么在这个文件里会重复出现,但是起始和终止位置稍有不同,可能因为编码出来的蛋白不同吧,所以同一基因也会被记录多次,所以存在部分同一位点交叉基因是同一基因的情况。
cat test1.txt | sort |uniq >test2.txt
当然也可以查看没有去除重复是有多少位点cat test.txt |awk -F "\t" '{sum+=$1-$2}END{print sum}'






回复 支持 反对

使用道具 举报

0

主题

4

帖子

52

积分

注册会员

Rank: 2

积分
52
发表于 2017-1-4 22:03:07 | 显示全部楼层
btrainee 发表于 2017-1-4 20:32
048-perl+R+python
一开始写的版本有点问题,把所有重复的坐标全去除了,没有考虑到不同的染色体。重新写 ...

您运行的结果和视频是一致的呀
回复 支持 反对

使用道具 举报

0

主题

12

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2017-1-4 22:23:39 | 显示全部楼层
003-python 发表于 2017-1-4 22:03
您运行的结果和视频是一致的呀

是我自己使用的同一个文件,但是用的两个不同的perl程序和那个python程序,3个程序运行出来的结果不同
回复 支持 反对

使用道具 举报

0

主题

4

帖子

52

积分

注册会员

Rank: 2

积分
52
发表于 2017-1-4 22:35:48 | 显示全部楼层
btrainee 发表于 2017-1-4 22:23
是我自己使用的同一个文件,但是用的两个不同的perl程序和那个python程序,3个程序运行出来的结果不同 ...

哦哦。。不好意思,只关注Python了
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-2-24 00:20 , Processed in 0.048668 second(s), 31 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.