搜索
楼主: Jimmy

生信编程直播第四题:多个同样的行列式文件合并起来

  [复制链接]

0

主题

4

帖子

50

积分

注册会员

Rank: 2

积分
50
发表于 2017-1-27 19:26:50 | 显示全部楼层
103-python-如梦
[Python] 纯文本查看 复制代码
def tj(lj):
    import os
    zd={}
    files=os.listdir(lj)        #提取所txt中的基因ID存入字典,这样就不会因不一致出问题
    for fl in files:
        f=open(lj+'/'+fl)
        for i in f:
            lb=i.split('\t')
            mh=lb[0].strip('\n')
            zd.setdefault(mh)
        f.close()
    zd2={}
    
    for fl in files:            #把每个文件中没有的基因添加在文件末尾并赋值为0,以便于遍历时规范
        zd2=zd.copy()
        f=open(lj+'/'+fl)       #这里需要注意,如果打开时后使用"a+"或"w+",那么就没法提取mh和去除字典中的相应条目
        for i in f:
            lb=i.split('\t')
            mh=lb[0].strip('\n')
            zd2.pop(mh)
        f.close()
        f=open(lj+'/'+fl,"a+")  #继续依次打开文件,把缺少的基因ID补全,为后面规范遍历做准备
        for i in zd2.keys():
            if not i=="":
                f.write(i+'\t'+'0'+'\n')
        f.close()
    
    for fl in files:            #最后一次遍历,把需要的数据存入字典
        f=open(lj+'/'+fl)
        for i in f:
            lb=i.split('\t')
            mh=lb[0].strip('\n')
            zd.update({mh:str(zd.get(mh)) +"★" +str(lb[1]).strip('\n')})
        f.close()

    for i in sorted(zd.keys(), reverse=True):
        print (str(i)+'\t'+str(zd.get(i).strip('None★')).replace("★",'\t'))

回复 支持 反对

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-1-28 20:59:22 | 显示全部楼层
写了一个TCGA版本的,先处理metadata
[Python] 纯文本查看 复制代码
import os
import sys
from collections import OrderedDict
import json
import re
import glob
 
 
gene_group=OrderedDict()
tcga_file=OrderedDict()
list_1=[]
normal=[]
tumor=[]
os.chdir("E:/z/files")
 
#处理TCGA metadata
with open('metadata.cart.2016-11-03T07_53_29.306026.json',"r") as f1:
    data=json.load(f1)
    for i in data:
        file_name=i["file_name"]
        file_name=re.sub("\.gz","",file_name)
        tcga_id=i["associated_entities"][0]["entity_submitter_id"]
        tcga_file[tcga_id]=file_name
        ID=tcga_id.split("-")
        if ID[3]=="11A":
            normal.append(tcga_id)              
        else:
            tumor.append(tcga_id)
list_1.extend(normal)
list_1.extend(tumor)
list_2=[]
for items in list_1:
    tcgafile=tcga_file[items]
    list_2.append(tcgafile)
##y=list_2[1]
#print(x,y)
def readfile(file):
    file_obj=open(file)
    try:
        while True:
            line=file_obj.readline().strip("\n")
            if not line:
                break
            arry=line.split("\t")
            if arry[0] in gene_group:
                gene_group[arry[0]].append(arry[1])
            else:
                gene_group[arry[0]]=[arry[1]]
    finally:
        file_obj.close()
    return gene_group
f3=open("TCGA-rnaseq-matrix.txt","w")
def main(): 
    for file in list_2:
        if glob.glob(file):           
            readfile(file)
    for gene_name,gene_expression in gene_group.items():
         
        f3.write(gene_name+"\t"+"\t".join(gene_expression)+'\n')
         
if __name__=="__main__":
    main()
回复 支持 反对

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-1-29 08:57:08 | 显示全部楼层
017中大普外 发表于 2017-1-28 20:59
写了一个TCGA版本的,先处理metadata[mw_shl_code=python,true]import os
import sys
from collections imp ...

赞,你懂得活学活用了,可以微信联系我,领取奖励!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-1-31 00:51:58 | 显示全部楼层
本帖最后由 017中大普外 于 2017-2-6 23:02 编辑
Jimmy 发表于 2017-1-29 08:57
赞,你懂得活学活用了,可以微信联系我,领取奖励!

谢楼主,差的还远,继续努力学习。。
回复 支持 反对

使用道具 举报

0

主题

6

帖子

53

积分

注册会员

Rank: 2

积分
53
发表于 2017-1-31 10:32:42 | 显示全部楼层
本帖最后由 snjiang12 于 2017-1-31 10:36 编辑

122-R-江
继续模仿学习中,讲讲自己体会。
先下载那个压缩文件,后面用git来合成一个,导入为tmp.txt。下面是R。看到楼主一步步地找,先以为是用melt函数,后面发现不是这个啊,原来是用cast,因为我们这个是数据框,所以用dcast函数。另外,对于这种大文件都是先取一小部分做测试,再整体使用的,任何函数循环也一样。下面是R中的代码:
[Python] 纯文本查看 复制代码
a<-read.table("tmp.txt",sep='\t',stringsAsFactors = F)
library(reshape2)
tmp<-a[1:20,]
names(tmp)<-tolower(names(tmp))#查起来是字符转向量的意思,如果没这个下面无法直接使用v1、v2
dcast(tmp,v2~v1)#v1是矩阵的名称,因而一个矩阵的所有样本都是这个名称,会自动合为一个的
#测试成功,下面用于整体
names(a)<-tolower(names(a))
result<-dcast(a,v2~v1)#结束


回复 支持 反对

使用道具 举报

0

主题

14

帖子

243

积分

中级会员

Rank: 3Rank: 3

积分
243
发表于 2017-2-1 17:11:03 | 显示全部楼层
本帖最后由 sleep在床上 于 2017-2-2 12:51 编辑

012-R+perl#sleep
仔细看了作业和答案,暂时留个包学习,一会在试试楼主的,让我接着用下DEseq试试后面的差异分析
[AppleScript] 纯文本查看 复制代码
my %total;
my $header = "Gene";
 
foreach my $i(0..@ARGV-1){
  $header.="\t$ARGV[$i]";
  open I,"$ARGV[$i]";
  while(<I>){
    chomp;
    my @tmp=split;
    $total{$tmp[0]}[$i]=$tmp[1];
  }
}
print $header,"\n";
foreach my $key(keys %total){
  foreach my $i(0..@ARGV-1){
    if(!$total{$key}[$i]){
      $total{$key}[$i]=0;
    }
  }
  print $key,"\t",join("\t",@{$total{$key}}),"\n";
}

以上用哈希存的时候title会被存,导致导入R的时候DEseq 报错会有nonnumberic
回复 支持 反对

使用道具 举报

13

主题

33

帖子

243

积分

版主

Rank: 7Rank: 7Rank: 7

积分
243
发表于 2017-2-1 22:22:33 | 显示全部楼层

092-R-python

今天学习了下,挺实用的,赞一个!
[AppleScript] 纯文本查看 复制代码
rm(list = ls())

## 首先在GSE48213_RAW目录里面生成tmp.txt文件
awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt

## awk 为linux 下批处理命令,输出当前文件名, $0表示整个记录--awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt

##去除表头---grep -v EnsEMBL_Gene_ID
##赋值给tmp.txt

## 然后把tmp.txt导入R语言里面用reshape2处理即可!
setwd('C:/Users/Administrator/Desktop/learn/lecture4/新建文件夹/tmp')
a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)

library(reshape2)
fpkm <- dcast(a,formula = V2~V1)
??dcast

row.names(fpkm)<-fpkm$V2
fpkm<-fpkm[,-1]

names(fpkm)<-unlist(lapply(names(fpkm),function(x){
  tmp<-strsplit(x,'_')[[1]][2]  
  tmp<-strsplit(tmp,'\\.')[[1]][1]
}))  ##体会[1]与[[1]]的区别

write.csv(fpkm,"56sample fpkm.csv")

回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-2-4 14:35:34 | 显示全部楼层
本帖最后由 旭日早升 于 2017-2-4 20:00 编辑

201-python-无名

本周作业是将每个样本的表达值合并为一个矩阵,提示注意的是可能gene排序顺序不同,我之前遇到过类似问题,应该是用R的merge做了一下,当时样本少做的也比较笨了,一个一个merge了,这次尝试用python搞定。

我的代码思路:首先初始化一个字典用于存入矩阵,依次读入每个文件的每一行,第一列是name,第二列是count,判断name是否在字典中,如果不在就加入字典并赋值count,如果在就直接把count加到对应元素上。遍历完后写入文件。问题,需要自己把文件名写入,不能自动化。于是根据东大神之前的指到,网上查找“python3 获得当前目录所有文件名”,使用os.listdir()可以解决。

我的代码:
[Python] 纯文本查看 复制代码
# Description: Cainiao4_v1
# Author: WKL
# Data: 2017-02-03

import re
import time 
import os
import gzip
from collections import OrderedDict
#calculate time
start = time.clock()

#change directory
path = "./GSE48213/"
os.chdir(path)

count_matrix = OrderedDict()

filenames = os.listdir()

for filename in filenames:
        with gzip.open(filename, 'rt') as fp:
                for line in fp:
                        genename = line.rstrip().split("\t")[0]
                        genecount = line.rstrip().split("\t")[1]
                        if genename in count_matrix:
                                count_matrix[genename] += ("\t" + str(genecount))
                        else:
                                count_matrix[genename] = str(genecount) 

with open("count_matrix.txt", 'w') as fp_out:
        for genename, genecount in count_matrix.items():
                fp_out.write("\t".join([genename, genecount]) + "\n")

end = time.clock()
print("used %s s" % str(end - start))


运行结果:矩阵文件
EnsEMBL_Gene_ID SUM229PE UACC812 。。。
ENSG00000000003 39.4250972624 26.5792803593 。。。
ENSG00000001167 111.1083686613 138.1403109667 。。。
。。。。。。




视频学习:学习了老师的讲解,老师更加注重结构和代码的可重复利用,函数的定义使得问题简单重复化,更加容易理解,虽然通过main调用还不是很理解,但是会继续学习。(补充,if __name__ == __main__是为了应对代码作为模块被导入时不需要执行主函数的情况。)附上结合老师的代码做的修改。
修改代码:
[Python] 纯文本查看 复制代码
import glob
import time 
import os
import gzip
from collections import OrderedDict
#calculate time
start = time.clock()

#change directory
path = "./GSE48213/"
os.chdir(path)

count_matrix = OrderedDict()

def Readfile(File):
        file_obj = gzip.open(File, 'rt')
        try:
                while True:
                        line = file_obj.readline().strip("\n")
                        if not line:
                                break
                        array = line.split("\t")

                        if array[0] in count_matrix:
                                count_matrix[array[0]].append(array[1])
                        else:
                                count_matrix[array[0]] = [array[1]]
        finally:
                file_obj.close()
        return count_matrix

def main():
        list_dirs = glob.glob("./*.txt.gz")
        for i in list_dirs:
                Readfile(i)
        with open("count_matrix_teacher.txt", 'w') as fp_out:
                for gene_name in count_matrix:
                        fp_out.write("%s\t%s\n" %(gene_name, "\t".join(count_matrix[gene_name])))

#作为模块使用代码时,该段代码将不执行
if __name__ == '__main__': 
        main()

end = time.clock()
print("used %s s" % str(end - start))


最后,本次作业是一个很好的实例应用,我们经常会遇到类似问题,这种训练是我们认识到编程的简便性。

回复 支持 反对

使用道具 举报

0

主题

6

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2017-2-4 18:36:38 | 显示全部楼层
python-60

今天下午补做了第3题和第4题,都是在未看老师的课程的情况下写的脚本,都成功了,还是有进步的,但是里面肯定有很多的代码可以优化,第4题的代码黏贴入下所示。

开始的时候我将数据导入了一个普通的字典中,输出之后发现首行到文件的中间去了,于是换了OrderedDict模块执行,所以行的顺序和原始文件的顺序会不一致,不知道这个在老师的课程里会不会得到解决。待我认真看了老师的课程之后再来改进代码。

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

exp_dic = collections.OrderedDict()
exps = []

path = '/Users/sunshine/practice/practice4-combine_files/'

f1 = open(path+'GSE48213_RAW'+'/'+'GSM1172844_184A1.txt','r')
lines = f1.readlines()
for line in lines:
    exp_id = line.strip('\n').split('\t')[0]
    exp_dic[exp_id] = []
f1.close()

f_output = open(path+'combined.txt','w')

F = os.listdir(path+'GSE48213_RAW')

for filename in F:
    if filename.startswith('.'):
        pass
    else:
        f = open(path+'GSE48213_RAW'+'/'+filename,'r')
        lines = f.readlines()
        for line in lines:
            line = line.strip('\n').split('\t')
            exp_id = line[0]
            exp = line[1]
            exp_dic[exp_id].append(exp)
        f.close()
    
for exp_id, exp in exp_dic.items():
    exp = '\t'.join(exp)
    f_output.write(exp_id+'\t'+exp+'\n')

f_output.close()
回复 支持 反对

使用道具 举报

1

主题

17

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
发表于 2017-2-5 21:59:11 | 显示全部楼层
Python-048
题目是合并表达量count文件,一般合并的count文件都是有相同的基因id,不会出现哪个文件中少一个gene的情况。但是,如果合并vcf文件,不同基因不同位点会发生不同的突变,所以python写的时候先遍历了一下所有文件的geneID存入字典,这样再进行合并的时候没有出现的gene相应位置直接步“-”好了。
代码如下:
[Python] 纯文本查看 复制代码
#! /usr/bin/env python3
# -*- coding:utf-8 -*-

import glob
from collections import OrderedDict

list_gene = OrderedDict()
text_filenames = glob.glob('~/home/btrainee/4_merge_count/test/*.txt') 

for filename in text_filenames:
#	print(filename)
	with open(filename,'r') as fp:
		for line in fp:
			lines = line.strip().split('\t')
			if not lines[0] in list_gene:
				list_gene[lines[0]] = ''

for filename in text_filenames:
	with open(filename,'r') as fp:
		for line in fp:
			lines = line.strip().split('\t')
			if not lines[0] in list_gene:
				list_gene[lines[0]] = list_gene[lines[0]] + '-\t'
			else:
				list_gene[lines[0]] = list_gene[lines[0]] + str(lines[1]) + '\t'
with open('result.xls','w') as fp_out:
	for key,value in list_gene.items():
		fp_out.write('\t'.join([key,value]) + "\n")

部分结果如下:
EnsEMBL_Gene_ID SUM149PT        CAMA1   HCC1599 ZR751   EFM192C MDAMB175VII     SUM225CWN       MX1     HCC3153 SUM229PE
ENSG00000000003 34.3229185245   21.1598896533   127.9318903236  31.8129987626   1.1394799588    44.5767736722   30.797926
ENSG00000001167 77.3913140961   381.1349474963  178.1828677921  68.4580223909   80.8550744211   122.2383747561  101.40240
ENSG00000005471 0.0000000000    0.8801682893    0.0000000000    0.5488649648    1.9199999586    0.0000000000    1.0061375
ENSG00000066629 0.0000000000    4.3341338195    4.2797094988    25.4937998700   114.0446564996  11.6454253395   22.957617
ENSG00000154258 0.0000000000    0.0312396434    0.0000000000    0.0000000000    0.0681193364    0.0000000000    0.0000000
ENSG00000154262 0.0000000000    0.0000000000    0.0000000000    0.0000000000    0.1215497979    0.0000000000    0.0000000
ENSG00000154263 0.0618638239    0.7370814775    0.2190705786    0.4767941843    0.0321547966    0.1207449604    0.1551086

回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.