How To Merge Cuffdiff Results To Make A Big Table?

Question: How To Merge Cuffdiff Results To Make A Big Table?

6.3 years ago by

United States

newDNASeqer • 670 wrote:

I started from 40 fastq files for running the RNASeq Pipeline (Tophat - cufflinks - cuffmerge), then I divided the 40 tophat folders/files to 8 processes of running cuffdiff (each process has 5 Tophat input files, accepted_hits.bam) with the same merge.gtf (from cuffmerge). I now need to put the cuffdiff results in a big table so as to compare the fkpm values among the 40 samples. The name of my tophat files are s1, s2, s3, ..., s39, and s40. Could someone tell me how to merge the cuffdiff results? thanks

cuffdiff • 2.9k views

ADD COMMENT • link •

modified 6.3 years ago by seidel ♦ 6.9k • written 6.3 years ago by newDNASeqer • 670

Which of the cuffdiff output files do you want to merge?

ADD REPLY • link written 6.3 years ago by Devon Ryan ♦ 93k

the file named gene_exp.diff

ADD REPLY • link written 6.3 years ago by newDNASeqer • 670

They should all have the same number of lines, so you can trivially script something in python/perl/whatever (I'm presuming that you just want the coordinates, gene name, and XLOC ID). BTW, is there a reason you're not just running one instance of cuffdiff with all of the files?

ADD REPLY • link written 6.3 years ago by Devon Ryan ♦ 93k

6.3 years ago by

seidel ♦ 6.9k

United States

seidel ♦ 6.9k wrote:

You could also aggregate your cuffdiff results fairly easily with R if you're familiar with it. Cuffdiff creates an output directory, so assuming directory names based on your tophat file names, the following would work to open files in 40 directories, grab out the column of interest, combine it to one table, and create an output file:

# declare a variable to hold fpkms
fpkm <- c()
# create directory names
dirnames <- paste("s", 1:40, sep=""))
# loop through all directories and grab fpkm columns
for( i in 1:length(dirnames) ){
  fname <- paste(dirnames[i], "/genes.fpkm_tracking",sep="")
  x <- read.table(file=fname, sep="\t", header=T, as.is=T)
  fpkm <- cbind(fpkm, x[,"FPKM"])
}
# name the columns
colnames(fpkm) <- dirnames
# name the rows, they're all in the same order
rownames(fpkm) <- x[,1]

write.table(fpkm, file="fpkm.txt", sep="\t", col.names=NA)

Purists would say you don't need for loops in R, but the code would be more obscure.

ADD COMMENT • link written 6.3 years ago by seidel ♦ 6.9k