SEQanswers (
-   Bioinformatics (
-   -   DE Seq (

AdrianJ217 08-30-2012 07:27 AM

DE Seq
Hi, I am new to using DE Seq. I have 2 files, output from htseq, containing the name of the genes with the # of reads next to them. One file has reads before induction of rna interference (uninduced) and another after (induced). So I want to compare the read count of each gene in both files and I believe I can do that with DE-Seq. When I read the manual, it only mentions inputting one file name in the command line in R. How do I input both files? This is my first time using DE Seq. Thank you1

dpryan 08-30-2012 10:09 AM

Just merge the two files. The easiest way for you to do that would be to open both in Excel and just copy and paste the second column from one of the files into the empty third column of the other. If you end up doing this sort of thing often then you can either write a small program to do that sort of thing or just do it in R.

Don't forget to remove the "no_feature..." and similar lines from the end.

AdrianJ217 08-30-2012 11:24 AM

Thank you
Great, thank you so much, and then can you tell me what commands I need to type into R in order to run it using my file?

dpryan 08-30-2012 12:31 PM

Save it as a CSV file and then follow the DESeq vignette, which is very well written.

AdrianJ217 08-31-2012 04:37 AM

dpryan, thank you so much. Just one more question. You said to copy the 2nd column of one to the empty 3rd column of the other, but why wouldnt I copy both columns of one to the other, to have 4 columns. Won't they mismatch if I just copy and paste the 2nd column?

ETHANol 08-31-2012 05:31 AM

Or you could run this R script to merge your HTSeq-count output files. Probably a little quicker then Excel if you have a lot of files. Anybody know how to do this with a BASH script??


#Ethan Ford
#for use with
#Put all your output files from HTseq-count in a directory
#change directories to that directory in BASH.
#Run script from dirctory with your files by typing in the command line (in BASH not in R) $ Rscript path/to/this/scirpt/multimergeRNA.R
#outputs a file called countstable.txt in the directory you ran the script from
#The column with your gene names must have the header "gene"

#To modify for other uses change:
#        1)'by="gene"' to the name of the column used as a reference for merging.
#        2)"countstable.txt" to the name of your desired outputfile

multimerge = function(mypath){
filenames=list.files(path=mypath, full.names=TRUE)
datalist = lapply(filenames, function(x){read.delim(file=x,header=T)})
Reduce(function(x,y) {merge(x,y,all=TRUE,by="gene",sort=FALSE)}, datalist)
merged <- multimerge(WD)
write.table(merged, "countstable.txt", sep="\t", quote=FALSE, row.names=FALSE, na="0")

dpryan 09-02-2012 01:46 AM

@AdrianJ217: They should be sorted the same to begin with. The DESeq vignette expects a file where the first column is the gene ID (or whatever the counts were done by) and subsequent columns are other samples. So, if you just copied both columns in then things would get more complicated.

@ETHANol: I don't have a bash script, but I did write a program in C to do that (it also changes the label on each column so it's easier to remember what is what). I'm sure I can send you a copy if you'd like, it makes making processing pipelines a bit easier.

AdrianJ217 09-02-2012 01:48 AM

dpryan, if you could send me a copy of that I'd really appreciate it. My email is:

Thank you so much!

ETHANol 09-02-2012 02:57 AM

dpryan, thanks!!! but in that I don't know any C at all, I don't think I could work it into my data analysis pipeline.

kopi-o 09-02-2012 03:41 AM


Anybody know how to do this with a BASH script??
You should be able to just loop over all files and use the 'join' command.

This of course assumes that the genes are always in the same order (which they are from HTSeq) but you may get into trouble if you have, e.g., HTSeq output using different versions of GTF files (like different ENSEMBL releases), when you will lose some genes from one of the files. That's why I think it's still better to join in R, Python etc where it is easier to check that nothing gets lost.

All times are GMT -8. The time now is 01:33 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.