搜索
查看: 689|回复: 18

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

[复制链接]

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1715
发表于 2016-12-29 09:19:19 | 显示全部楼层 |阅读模式
相信用过htseq-count的朋友都知道,它是分开对每个样本计算所有的基因表达量,所以会生成一个个独立的文件,我用perl脚本模仿它的结果如下:
$ head a.txt
gene_1  178
gene_2  692
gene_3  486
gene_4  666
gene_5  395
gene_6  48
gene_7  926
gene_8  733
gene_9  660
gene_10 578

第一列是基因,第二列是该基因的counts值,共有a~z这26个样本的counts文件,需要合并成一个大的行列式,这样才能导入到R里面做差异分析,如果手工用excel表格做,当然是可以的,但是太麻烦,如果有500个样本,正常人都不会去手工做了,需要编程。

生成测试文件的代码如下:
[AppleScript] 纯文本查看 复制代码
#首先新建文件tmp.sh 输入这个代码:
perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'
## 然后用perl脚本调用这个tmp.sh文件:
perl -e 'system(" bash tmp.sh >$_.txt") foreach a..z'
##这样就生成了a~z这26个样本的counts文件


用shell或者perl或者python,设置R语言都可以做,但是各有优缺点,而且如果每个样本的基因顺序并不一致,这时候你应该怎么做呢?

实际需求如下: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的不同,分组做差异分析。
paper是:https://www.ncbi.nlm.nih.gov/pubmed/24176112

输出的表达矩阵,如下所示:


请务必做出此题,很重要!


先给一下shell结合R语言的做法:
[Python] 纯文本查看 复制代码
## 首先在GSE48213_RAW目录里面生成tmp.txt文件
awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt
## 然后把tmp.txt导入R语言里面用reshape2处理即可!
setwd('tmp/GSE48213_RAW/')
a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)



本帖子中包含更多资源

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

x



上一篇:看了一篇HATs和HDACs的文章
下一篇:结构变异SV是什么?
回复

使用道具 举报

1

主题

9

帖子

104

积分

注册会员

Rank: 2

积分
104
发表于 2017-1-8 20:13:57 | 显示全部楼层
本帖最后由 disheng 于 2017-1-8 20:37 编辑

043-perl-R-shell-涤生

学习了二楼yydBIG的二维哈希,真的太棒了。参照他的思想,用哈希也写了个
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl

my %hash;
my $header = "Gene";
foreach my $i ( 0..@ARGV-1){
        $header .= "\t$ARGV[$i]";
        open I,"$ARGV[$i]";
        while(<I>){
                chomp;
                my($gene_num,$num) = split;
                $hash{$gene_num} .= "$num\t";
        }
}
print $header,"\n";
foreach $i (keys %hash){
        print "$i\t$hash{$i}\n";
}


本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

1

主题

2

帖子

21

积分

新手上路

Rank: 1

积分
21
发表于 2016-12-29 11:12:00 | 显示全部楼层
本帖最后由 yydBIG 于 2016-12-29 11:13 编辑

贴一个perl实现的版本,不用排序,没出现的用0代替。使用方法:
perl zuhe.pl file1 file2 ...
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
#2010-12-24
#yangyd

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";
}


回复 支持 1 反对 0

使用道具 举报

0

主题

6

帖子

69

积分

注册会员

Rank: 2

积分
69
发表于 2016-12-29 09:32:23 | 显示全部楼层
大多数情况下基因顺序是一致的
回复 支持 反对

使用道具 举报

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1715
 楼主| 发表于 2016-12-29 09:41:44 | 显示全部楼层
liyan 发表于 2016-12-29 09:32
大多数情况下基因顺序是一致的

这种情况下是一致的,但是,如果呢,比如其它软件,其它情况下,如果你需要合并vcf文件,坐标就不一致了,怎么办,能hold住吗?
回复 支持 反对

使用道具 举报

0

主题

1

帖子

18

积分

新手上路

Rank: 1

积分
18
发表于 2016-12-29 10:26:34 | 显示全部楼层
本帖最后由 l0o0 于 2016-12-29 10:49 编辑

这里我觉得用R语言里面的merge来做感觉会比perl或python建字典方便许多。
这里不方便贴代码,可以看下github上的链接https://github.com/l0o0/bio-analysis-kit/blob/master/merge.R,默认文件第一行是表头。
[AppleScript] 纯文本查看 复制代码
#!/bin/env Rscript

args = commandArgs(T)

if (length(args) != 4){
    print('USAGE: Rscript merge.R pattern merge_colname input_dir output_file')
    quit()
}


files = list.files(args[3],pattern = paste(args[1],'$',sep=''))

a <- files[1]
atable <- read.table(a, header=T, check.names=F, sep='\t')
samplea = strsplit(basename(a), '.', fixed=T)[[1]][1]
colnames(atable) = c(args[2], samplea)

for (i in files[-1]){
	print (i)
	sampleb = strsplit(basename(i), '.', fixed=T)[[1]][1]
    btable <- read.table(i, header=T, check.names=F, sep='\t')
    colnames(btable) = c(args[2], sampleb)
	print(colnames(btable))
	#names(btable) <- c('gene_id', fpkm)
	atable <- merge(atable, btable, by = args[2], all =TRUE)
}

row.names(atable) = atable[,1]
atable = atable[,-1]
print(colnames(atable))
write.table(atable, sep='\t', file=args[4], col.names=NA,quote=FALSE, na='NA')

Rscript merge.R file-pattern  geneid $PWD test.txt

回复 支持 反对

使用道具 举报

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1715
 楼主| 发表于 2016-12-29 10:31:33 | 显示全部楼层
l0o0 发表于 2016-12-29 10:26
这里我觉得用R语言里面的merge来做感觉会比perl或python建字典方便许多。[mw_shl_code=applescript,true]ar ...

我给你示范了一下如何排版,你继续修改一下你的代码,稍微正常一点,编辑器里面有插入代码选项
回复 支持 反对

使用道具 举报

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1715
 楼主| 发表于 2016-12-29 15:46:05 | 显示全部楼层
yydBIG 发表于 2016-12-29 11:12
贴一个perl实现的版本,不用排序,没出现的用0代替。使用方法:
perl zuhe.pl file1 file2 ...
[mw_shl_cod ...

你的代码很简练,赞,虽然我没有测试过,但是跟我的思路是一样的
回复 支持 反对

使用道具 举报

3

主题

7

帖子

41

积分

新手上路

Rank: 1

积分
41
发表于 2017-1-4 22:44:56 | 显示全部楼层
htseq-count出来一样的行的,我会用下面这个命令来做
[Shell] 纯文本查看 复制代码
paste file1 file2 file3 | awk '{printf $1"\t";for(i=2;i<=NF;i=i+2)printf $i"\t";print $i}'
回复 支持 反对

使用道具 举报

274

主题

473

帖子

1715

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1715
 楼主| 发表于 2017-1-4 23:33:07 | 显示全部楼层
猪兔子2号 发表于 2017-1-4 22:44
htseq-count出来一样的行的,我会用下面这个命令来做
[mw_shl_code=shell,true]paste file1 file2 file3 |  ...

你这个也很赞
回复 支持 反对

使用道具 举报

1

主题

38

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-1-6 10:54:28 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-6 11:12 编辑

059pyth+R 冷洋
意思是要把所有的文本文件的内容合成一个文本文件吗?
[Python] 纯文本查看 复制代码
#coding=utf-8
import os
import os.path
g=open('D:/one.txt','w')
path='D:/56-breast-cancer-celline-RNA-seq'
files=os.listdir(path)#得到path路径下的所有文件名称
for file in files:
    f=open(path+'/'+file)#逐个打开文件
    for i in f:
        g.write(i+'\n')#将每个文件的内容写入one.txt
g.close()

回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2017-1-22 11:49 , Processed in 0.242632 second(s), 33 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.