搜索
楼主: Jimmy

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

  [复制链接]

0

主题

9

帖子

79

积分

注册会员

Rank: 2

积分
79
发表于 2017-2-6 21:20:48 | 显示全部楼层
PYTHON-010
[AppleScript] 纯文本查看 复制代码
import glob
import os
os.chdir('C:/Users/Administrator/Desktop/')

d={}
#先定义一个字典,这个字典将文件的内容都存进去,格式是,key:基因名字,value:值的数组

def readfile(file):
    file_obj=open (file) 
 #先把文件打开 
    try:
        while True:
            line=file_obj.readline().strip('\n')
#上面两行代码,一行一行读取文件,每一行一个循环
            if not line:
                break
#如果是空的话,就跳出这循环直接跳到finally,结束这个文件的读取
            array=line.split()
            if array[0] in d:
                d[array[0]].append(array[1])
#字典里面如果没有这个gene名字,后面加上这个gene的值
            else:
                d[array[0]]=[array[1]]
#字典里如果没有,则将gene名字加进去,然后定义值,注意值是列表的形式
    finally:
            file_obj.close()
    return d
#最后返回字典值


def main():
    list_dirs=glob.glob('./*.txt')
#读取当先文件夹里的所有的txt文件
    for i in list_dirs:
        readfile(i)
#使用函数一个一个的读入字典中
    for gene_name in d:
        print(gene_name,'\t'.join(d[gene_name]))
# join() 方法将值一个一个的用tab键分隔开
        

        
if __name__=='__main__':
    main()
#最后一行不大懂
回复 支持 反对

使用道具 举报

0

主题

4

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-2-6 21:37:10 | 显示全部楼层
本帖最后由 zhongzheng 于 2017-2-6 21:38 编辑

212-perl-zhongzheng
两个思路都试了一下,下面代码是哈希中嵌套数组变量:
[Perl] 纯文本查看 复制代码
#! /usr/bin/perl
use strict;
use warnings;

my %total;
my $header = "Gene";
foreach my $i (0..@ARGV-1){
    $ARGV[$i] =~ /(.*?)_(.*?).txt/;
    $header .= "\t$2";
  open I, "$ARGV[$i]";
  while (<I>){
    chomp;
        my @tmp = split;        
        $total{$tmp[0]}[$i] = $tmp[1];
  }
}

print $header, "\n";
foreach my $key (sort keys %total){
    foreach my $i (0..@ARGV-1){
          if (!$total{$key}[$i]){
            $total{$key}[$i]=0;
          }
        }        
        print $key,"\t",join ("\t",@{$total{$key}}),"\n";
}




本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-2-6 22:10:25 | 显示全部楼层
074-python
用的是论坛帖子中做WGCNA分析的56个sample,主要是理解利用字典存储结果,且将数组作为字典的值这个过程。
[Python] 纯文本查看 复制代码
#!/usr/bin/python
#-*-coding:utf-8-*-

import glob  #这个包可以查找符合特定规则的文件路径名

Mydic = {} #定义一个储存结果的字典

def Readfile(File):  #定义一个函数
		file_obj = open(File)   #File这里是个形参
		try:                 #try..finally 是python读取文件的常用格式
		    while True:
		    	  line = file_obj.readline().strip("\n") #readline返回一个list,每行是一个list,strip去掉每行后的换行符
		    	  if not line:
		    	  	  break
		    	  array = line.split("\t")  #每一行按照\t分割成一个数组
		    	  
		    	  if array[0] in Mydic :
		    	  	  Mydic[array[0]].append(array[1])  #键对应一个字符串,值对应一个数组。当字典里存在基因名的键时,把该基因的表达量加到该键对应的值中
		    	  else:
		    	      Mydic[array[0]] = [array[1]] #如果genename(array[0])第一次出现,把基因的表达量作为数组的第一个元素
		finally:
			  file_obj.close()
		return Mydic
		
def main():
	  list_dirs = glob.glob("./*.txt")  #glob.glob返回所有匹配的文件路径列表
	  for i in list_dirs:
	  	  Readfile(i)
	  for gene_name in Mydic:
	  	  print "%s\t%s" % (gene_name,"\t".join(Mydic[gene_name])) #用圆括号元组里的值替换前面的%s
	  	  
if __name__ == '__main__':
	  main()

回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-2-7 11:18:38 | 显示全部楼层
end2end 发表于 2017-1-6 19:06
I think pandas can make our job easier!

[mw_shl_code=python,true]import pandas

你好!学习你的代码运行之后,最后生成的csv文件的样本名是原样本名分割之后的,请教一下哪一行代码是对文件名进行处理?
回复 支持 反对

使用道具 举报

0

主题

13

帖子

93

积分

注册会员

Rank: 2

积分
93
发表于 2017-2-7 22:08:39 | 显示全部楼层
本帖最后由 tingting 于 2017-2-7 22:15 编辑

121-R&python-婷

第四题又有一点小进步,就是完全没看视频做了题目,虽然写得非常差,但是勉强得到了结果,做完后再看老师视频进行优化改进。

第四题题目:多个同样的行列式文件合并起来(合并不同样品测序所得的基因FPKM值)
先用下面的脚本检测一下字典是否有序:
gene_id = {}      #后面又把这里换成gene_id = OrderedDict()
gene_id['gene_1'] =258
gene_id['gene_4'] =238
gene_id['gene_2'] =268
gene_id['gene_3'] =248
print gene_id
lst = sorted(gene_id)
print lst

结果发现:sorted以后是这样的['gene_1', 'gene_2', 'gene_3', 'gene_4']
而{}字典的排序本来就是这样{'gene_1': 258, 'gene_2': 268, 'gene_3': 248, 'gene_4': 238},
OrderedDict却是这样的([('gene_1', 258), ('gene_4', 238), ('gene_2', 268), ('gene_3', 248)]),所以本题需要用普通字典{}
最后,发现上面这个是白做的,不需要判断字典的有序性。不过通过这个小测试明白了普通字典与OrderedDict的区别。


下面是最开始写的程序,虽然很差劲,但是想留作纪念,希望以后越来越好!

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

args = sys.argv
gene_dict = {}
#print os.listdir('F:test_data/training_tree_4/')
dir_files = os.listdir('F:test_data/training_tree_4/')           #这里不知道是不是可以替换成args[1]?
#print dir_files
os.chdir('F:test_data/training_tree_4/')                         #这里不知道是不是可以替换成args[1]?
for i in dir_files:
    with open(i) as f:
        for line in f:
            line = line.rstrip()
            lst = line.split('\t')
            if lst[0] not in gene_dict.keys():
                gene_dict[lst[0]]= [int(lst[1])] 
            else:
                gene_dict[lst[0]].append(int(lst[1]))
for gene_id, gene_lst in gene_dict.items():
    print gene_id, gene_lst


下面是结果:
gene_1 [224, 23, 682, 131, 105, 41, 390, 794, 567, 597, 898, 66, 8, 741, 244, 715, 241, 499, 628, 125, 186, 950, 959, 259, 981, 524]
gene_2 [225, 722, 931, 185, 736, 194, 441, 172, 502, 222, 913, 635, 282, 101, 798, 259, 963, 96, 694, 233, 319, 160, 421, 467, 902, 862]
gene_3 [584, 683, 898, 719, 976, 952, 566, 674, 47, 881, 690, 556, 463, 456, 336, 282, 940, 823, 686, 172, 851, 975, 692, 288, 194, 291]
gene_4 [41, 954, 230, 867, 603, 43, 430, 669, 137, 735, 617, 821, 133, 669, 696, 116, 615, 438, 485, 504, 767, 94, 301, 170, 507, 373]
gene_5 [48, 89, 837, 88, 474, 613, 780, 194, 743, 235, 100, 266, 666, 412, 852, 177, 798, 443, 916, 325, 402, 815, 350, 592, 501, 904]
gene_6 [657, 156, 16, 781, 292, 854, 411, 685, 228, 565, 88, 474, 848, 566, 646, 348, 408, 50, 948, 785, 513, 628, 540, 564, 496, 87]
gene_7 [883, 162, 588, 442, 615, 213, 935, 807, 666, 507, 835, 520, 818, 741, 503, 177, 193, 205, 374, 889, 453, 40, 364, 813, 202, 647]
gene_8 [638, 862, 605, 159, 815, 302, 38, 424, 15, 649, 309, 819, 504, 972, 107, 205, 481, 721, 660, 313, 865, 946, 423, 990, 799, 80]
gene_9 [534, 335, 507, 705, 821, 137, 930, 150, 966, 259, 291, 620, 502, 194, 699, 323, 952, 155, 262, 401, 817, 391, 197, 825, 880, 536]

不知道为什么跑这个脚本第二遍的时候结果就出不来了,出来的是错误(错误信息如下),只能关掉jupyter重启就又可以跑出一次结果,不知道哪里有问题,如果哪位老师看到了帮忙解答一下,谢谢!

---------------------------------------------------------------------------
WindowsError                              Traceback (most recent call last)
<ipython-input-2-be15552a607f> in <module>()     21 gene_dict = {}     22 #print os.listdir('F:test_data/training_tree_4/')---> 23 dir_files = os.listdir('F:test_data/training_tree_4/')     24 #print dir_files     25 os.chdir('F:test_data/training_tree_4/')WindowsError: [Error 3] : 'F:test_data/training_tree_4/*.*'



回复 支持 反对

使用道具 举报

2

主题

34

帖子

777

积分

高级会员

Rank: 4

积分
777
发表于 2017-2-9 14:22:21 | 显示全部楼层
本帖最后由 x2yline 于 2017-9-13 19:45 编辑

077-R-python-shell x2yline

1.数据下载
chrome浏览器进入网址https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213 后点击底部的http超链接进行下载,下载好的文件GSE48213_RAW.tar直接放在G盘下,不用进行解压

2.python实现
先贴上代码:
[Python] 纯文本查看 复制代码
import os

def untar(fname, untar_path):
    ''' 函数输入为fname为tar或gz文件的绝对路,
    untar_path为解文件放置路径,不存在时可以自动新建文件夹,
    utar_path若为空,则默认为fname的同级新建目录'''
    import tarfile
    if not untar_path:
        untar_path = ''.join(fname.split('.')[:-1])
    t = tarfile.open(fname, 'r')
    t.extractall(path = untar_path)
    t.close()

def merge_array(file_list, dirs = '.', na_fill='' ):
    '''函数输入file_list为要和并的距阵文件名列表,可以是txt或gz文件
    dirs为文所在的路径,默认为当前路径
    输出为有序字典,keys为行名,values为对应行的数值列表'''
    my_dict = {}
    counts = -1
    for i in file_list:
        print(i)
        counts += 1
        if i[-3:] == '.gz':
            import gzip
            with gzip.open(os.path.join(dirs, i),'r') as f:
                for line in f:
                    item = line.decode("utf-8").split('\t')
                    if len(item) <=1 :
                    	print(item)
                    	break
                    if item[0] in my_dict:
                        if (counts - len(my_dict[item[0]])):
                            print("Warning: a record missing in previous files: " + item[0] + "\n")
                        my_dict[item[0]] += (counts - len(my_dict[item[0]]))*[na_fill,] + [item[1].strip()]
                    else:
                        my_dict[item[0]] = counts*[na_fill,] + [item[1].strip()]
                        if counts:
                            print("Warning: a record missing in previous files: " + item[0] + "\n")
        else:
            with open(os.path.join(dirs, i),'r') as f:
                for line in f:
                    item = line.split('\t')
                    if item[0] in my_dict:
                        if (counts - len(my_dict[item[0]])):
                            print("Warning: a record missing in previous files: " + item[0] + "\n")
                        my_dict[item[0]] += (counts - len(my_dict[item[0]]))*[na_fill,] + [item[1].strip()]
                    else:
                        my_dict[item[0]] = counts*[na_fill,] + [item[1].strip()]
                        if counts:
                            print("Warning: a record missing in previous files: " + item[0] + "\n")
    return my_dict

def write_dict_to_csv(my_dict, file_list, path = '.\\merged.csv'):
    with open(path,'w') as f:
        # 写上表头
        f.write(','+','.join([j.split('.')[0] for j in file_list]) + '\n')
        for i in my_dict.keys():
            content = i + ','
            content += ','.join(my_dict[i]) + '\n'
            f.write(content)
    return 0

# fname为GSE48213_RAW.tar所在绝对路径
# untar_path为解压后文件所在路径,可任意设定
fname = r'G:\GSE48213_RAW.tar'
untar_path = r'G:\GSE48213_RAW'

untar(fname, untar_path)
file_list = []
for j in os.listdir(r'G:\GSE48213_RAW'):
    if j[-3:] == '.gz':
        file_list.append(j)
        
my_dict = merge_array(file_list, dirs = untar_path)
content = write_dict_to_csv(my_dict, file_list=file_list, path = 'E:\\merged_array.csv')


最后结果E:\\merged_array.csv'截图为:

相关说明:
(1)OrderedDict是有序的字典,它的顺序是根据键和值创建的时间进行排序的,所以第一个文件中Gene_ID的排序决定字典中Gene_ID的排序,但是由于文件中Gene_ID的排序并不是按照编号由小到大进行排序的,所以最终要实现Gene_ID的排序都要用到排序函数,这样OrderedDict就没有什么用了,因此用普通的字典和排序字典都需最后统一排序,我的代码中就只用了普通字典。
(2)代码中考虑了不同矩阵中行数不同的情况
(3)开始下载.tar文件时出现解压错误,最后换了一种方式(下载地址和下载器)下载数据而得以解决。

3.R语言实现

把GSE48213_RAW.tar解压到E盘的E:\\GSE48213_RAW目录下,生成56个.gz文件

运行以下代码即可
[AppleScript] 纯文本查看 复制代码
merge_file_by_column1 <- function(fname_list) {
  #函数输入为目标文件名列表,gz文件或txt文件
  #输出为按第一列合并的数据框
  pb <- txtProgressBar(1, length(fname_list)-1, style=3)
  merged_dat <- read.table(gzfile(fname_list[1]), stringsAsFactors = F, header = T)
  for (i in (1:(length(fname_list)-1))){
    setTxtProgressBar(pb, i)
    tmp_dat <- read.table(gzfile(fname_list[i+1]), stringsAsFactors = F, header = T)
    merged_dat <- merge(merged_dat, tmp_dat, by = colnames(tmp_dat)[1], all=T)
  }
  close(pb) 
  return(merged_dat)
}

# 先把GSE48213_RAW.tar解压为.gz文件(不用解压为TXT文件)
# 确保E:\\GSE48213_RAW目录下只含有目标文件(即只有56个gz文件或者只有56个txt文件)
setwd('E:\\GSE48213_RAW')

merged_dat <- merge_file_by_column1(fname_list=list.files())
View(merged_dat)

最后结果得到一个根据第一列合并的datafarme变量




回复 支持 反对

使用道具 举报

6

主题

23

帖子

201

积分

中级会员

Rank: 3Rank: 3

积分
201
发表于 2017-2-9 20:17:19 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-2-9 20:28 编辑

1. 首先下载数据
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213 里面有56个文件(ftp://ftp.ncbi.nlm.nih.gov/geo/s ... pl/GSE48213_RAW.tar),需要合并成一个表达矩阵,来根据cell-line的不同,分组做差异分析。
[AppleScript] 纯文本查看 复制代码
wget [url=ftp://ftp.ncbi.nlm.nih.gov/geo/s]ftp://ftp.ncbi.nlm.nih.gov/geo/s[/url] ... pl/GSE48213_RAW.tar
tar -zxvf GSE48213_RAW.tar
for i in GSM* ;do  gzip -d *;done
ls 


2. 介绍一下glob这个包
用于查找符合特定规则的文件路径名。有3个匹配符:
* : 匹配0个或多个字符
? : 匹配单个字符
[] :匹配指定范围内的字符

glob.glob
返回所有匹配文件路径的列表
例如:
[AppleScript] 纯文本查看 复制代码
>>> import glob
>>> print glob.glob("*.txt")


glob.iglob
获取一个可遍历的对象,可以逐个获得匹配的文件路径名
[AppleScript] 纯文本查看 复制代码
>>> f = glob.iglob("*.txt")
>>> for i in f:
...     print i



3. try/else语句
try语句的完整形式:
try:
    <···>
except <  >:
    <···>
else:
    <···>
except用于try代码发生异常时执行,else用于正常执行try后执行

4.try/except/finally语句
无论try语句中是否有异常,finally语句都会执行!
5.格式化输出
%格式符
s,获取传入对象的返回值,并将其格式化到指定位置
[AppleScript] 纯文本查看 复制代码
>>> s = 'aaa %'
>>> print(s)
aaa %
>>> s1 = 'aaa %s'%('bbb')
>>> print(s1)
aaa bbb
>>> 


6.需要合并成一个表达矩阵
[AppleScript] 纯文本查看 复制代码
#! /bin/python
import os
import glob
from collections import OrderedDict
MyDiction = OrderedDict()
def Readfile(File):
        file_obj = open(File)
        try:
                while True:
                        line = file_obj.readline().strip("\n")
                        if not line:
                                break
                        array = line.split("\t")
                        if array[0] in MyDiction:
                                MyDiction[array[0]].append(array[1])
                        else:
                                MyDiction[array[0]] = [array[1]]
        finally:
                file_obj.close()
        return MyDiction
def main():
        list_dirs = glob.glob("./*.txt")
        for i in list_dirs:
                Readfile(i)
        with open("count_matrix.txt", 'w') as f_out:
                for gene_name in MyDiction:
                        f_out.write("%s\t%s" %(gene_name,"\t".join(MyDiction[gene_name])))
if __name__ == '__main__':
                main()
        

生成:
回复 支持 反对

使用道具 举报

2

主题

9

帖子

235

积分

中级会员

Rank: 3Rank: 3

积分
235
发表于 2017-2-10 10:33:55 | 显示全部楼层
本帖最后由 wangfangce 于 2017-2-10 10:46 编辑

206-Perl-策
因为perl是一行一行读取所以先要让数据按相同基因的放在一起,再按样本顺序排列,所以在tmp.txt的基础上加上一列为基因名连接样品的列,再sort这一列。就可以依次输出了。
cat tmp.txt|perl -alne'print "$F[1]$F[0]\t$F[0]\t$F[1]\t$F[2]"'>2_tmp.txt
然后
sort  2_tmp.txt >sorted.txt
再用perl脚本处理如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
open IN,"sorted.txt";
$last="ENSG00000000003";
print $last;
while(<IN>){
@F=split; 
$hash{$F[1]}=1;
if ($F[2] eq $last){
print "\t$F[3]";
}else{
print "\n$F[2]\t$F[3]";
$last=$F[2];
}
}
print "\n\t";
foreach $key(sort keys %hash){
print "$key\t";
}
结果截图
样品信息在最后一行,在EXCEL里调整到第一行就行了。 好吧 我的方法太麻烦了 。 看看别人的学习学习,

本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

5

帖子

115

积分

注册会员

Rank: 2

积分
115
发表于 2017-2-10 19:11:22 | 显示全部楼层
看视频前
先采用测试数据,进一步简化到只是处理2个数据,a.txt, b.txt;
思路:
采用hash识别ab文件中gene (key)和count (value),利用sort keys 遍历,exists keys 函数,识别相同的gene(key)
问题,如果有26个文件,甚至更多,怎么办?一个个的写,还是很麻烦哎
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;


open IN_a, 'a.txt' or die $!;
open IN_b, 'b.txt' or die $!;

my ($gene, $count, %hash_a, %hash_b);

while (<IN_a>){
	chomp;
	($gene, $count)=split/\t/;
	$hash_a{$gene}=$count;
		}
	
while (<IN_b>){
	chomp;
	($gene, $count)=split/\t/;
	$hash_b{$gene}=$count;
	}
foreach (sort keys %hash_b){
    if (defined $hash_a{$_}){ 
    	# or exists $hash_a{$_}
	    print "$_\t$hash_a{$_}\t$hash_b{$_}\n";
}
}
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2017-2-10 23:34:24 | 显示全部楼层
本帖最后由 chapman 于 2017-5-20 12:36 编辑

[Python] 纯文本查看 复制代码
import glob
mydict = {}

def readfile(file):
    f = open(file)
    for lines in f:
        line = lines.rstrip().split('\t')
        gene_name = line[0]
        gene_info = line[1]
        if gene_name in mydict:
            mydict[line[0]].append(line[1])
        else:
            mydict[line[0]] = [line[1]]

    f.close()


def main():
    list_dir = glob.glob('./*.txt')
    for i in list_dir:
        readfile(i)
    ft = open('out_1.txt', 'w')
    for genename,info in mydict.items():

        ft.write(genename + '\t' + "\t".join(info) + '\n')

if __name__ == "__main__":
    main()

0、050-Python-
1、我发现这个算法好慢,再看看有没有算法好点的。(已经找到,用了for in 逐行读取,没有readline);
2、我发现我过个年,我又忘记好多了,要赶紧找回来。

回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-13 03:57 , Processed in 0.035885 second(s), 20 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.