SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Reply
 
Thread Tools
Old 08-30-2012, 07:27 AM   #1
AdrianJ217
Junior Member
 
Location: Israel

Join Date: Aug 2012
Posts: 9
Default 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
AdrianJ217 is offline   Reply With Quote
Old 08-30-2012, 10:09 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

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.
dpryan is offline   Reply With Quote
Old 08-30-2012, 11:24 AM   #3
AdrianJ217
Junior Member
 
Location: Israel

Join Date: Aug 2012
Posts: 9
Default 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?
AdrianJ217 is offline   Reply With Quote
Old 08-30-2012, 12:31 PM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Save it as a CSV file and then follow the DESeq vignette, which is very well written.
dpryan is offline   Reply With Quote
Old 08-31-2012, 04:37 AM   #5
AdrianJ217
Junior Member
 
Location: Israel

Join Date: Aug 2012
Posts: 9
Default

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?
AdrianJ217 is offline   Reply With Quote
Old 08-31-2012, 05:31 AM   #6
ETHANol
Senior Member
 
Location: Western Australia

Join Date: Feb 2010
Posts: 308
Default

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??

Code:
#!/usr/bin/R
#multimergRNA.R 
#Ethan Ford
#for use with ezDESeq.sh
#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")
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 09-02-2012, 01:46 AM   #7
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

@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.
dpryan is offline   Reply With Quote
Old 09-02-2012, 01:48 AM   #8
AdrianJ217
Junior Member
 
Location: Israel

Join Date: Aug 2012
Posts: 9
Default

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

Thank you so much!
AdrianJ217 is offline   Reply With Quote
Old 09-02-2012, 02:57 AM   #9
ETHANol
Senior Member
 
Location: Western Australia

Join Date: Feb 2010
Posts: 308
Default

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.
__________________
--------------
Ethan
ETHANol is offline   Reply With Quote
Old 09-02-2012, 03:41 AM   #10
kopi-o
Senior Member
 
Location: Stockholm, Sweden

Join Date: Feb 2008
Posts: 319
Default

Quote:
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.
kopi-o is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 09:03 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO