SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
To convert BAM file to BED file Anjali Bioinformatics 10 04-28-2014 05:43 AM
Gene 12 Column BED file obifro Bioinformatics 2 10-21-2013 10:09 AM
TopHat file.bam file.bed join Trudy Bioinformatics 1 05-21-2013 12:59 PM
find overlaps/common in multiple bed file epi Bioinformatics 11 02-05-2013 06:47 AM
what is the file size for a 30X human genome sequencing file, raw and BAM? RNA-seq Illumina/Solexa 2 04-15-2011 12:27 PM

Reply
 
Thread Tools
Old 02-25-2013, 12:46 PM   #1
swaraj
Member
 
Location: Naples, Italy

Join Date: Feb 2012
Posts: 50
Default Raw counts of 12 column bed file (against multiple BAM)

Hello Everyone,

I am sharing a perl script to get raw read counts from mulitple BAM files, for a given bed file. The script requires bedtools and samtools to be preinstalled in the system.

I hope it is useful for users. An alternative solution exists in python
(http://seqanswers.com/forums/showthread.php?t=18182)
Attached Files
File Type: pl bed_count.pl (5.7 KB, 19 views)
swaraj is offline   Reply With Quote
Old 09-23-2013, 07:51 AM   #2
jason_ARGONAUTE
Member
 
Location: china

Join Date: Aug 2013
Posts: 14
Default

thanks so much for ur script.
But you didn't give the usage info.

I hv a .gtf file download from EMBL for D.mel now, should i convert it to .bed in advance?
and the usage is?
perl bed_count.pl annotated.bed dir_of_my_bamfiles tmpcount/ my.exon my.count
jason_ARGONAUTE is offline   Reply With Quote
Old 09-23-2013, 08:30 AM   #3
swaraj
Member
 
Location: Naples, Italy

Join Date: Feb 2012
Posts: 50
Default

Dear Jason,

A much quicker solution is to

1. Convert GTF to BED. A nice script is found here http://ea-utils.googlecode.com/svn/t...lipper/gtf2bed

2. Convert the BED file to exon coordinates (BEDTools package).
bed12toBed6 -i filename.bed > filename.exons

3. Get read counts from multiple BAM files (BEDTools package)
multiBamCov -bams aln.1.bam aln.2.bam -bed filename.exons > filename.counts
(Remember to build index of each BAM file using samtools index, else multiBamCov wont work)

4. Load the count file in R environment.

filename <- read.delim("filename.counts", header=FALSE,sep = "\t", stringsAsFactors =FALSE)

5. Sum the exons into gene counts using apply function in R
input <- filename[,7:ncol(filename)]
names <- filename$V4

filename.gene.count<-apply(input,2,function(x){tapply(x,names,function(y){sum(y)})})

6. Write the count file as a text file
write.table(filename.gene.count, file = "filename.gene.count", quote = FALSE, row.names=FALSE, col.names=FALSE)

The order of the counts for samples will be the same order in which you give the BAM files. My previous script is easier to implement but little slower. If you open the script in a text editor you can give the name of the directory of BAM files and the BED file. Hope this helps.
swaraj is offline   Reply With Quote
Old 09-24-2013, 01:37 AM   #4
jason_ARGONAUTE
Member
 
Location: china

Join Date: Aug 2013
Posts: 14
Default

Dear swaraj,

It helps me a lot. Thanks again^^
But i hv encountered a now problem now when adopting the new solution. I hv got 23017 exons in total while 14797 genes are mapped to the genome by Tophat & Cufflinks. and 14797 is a proper count for fruit fly. I notice that ur solution use exon coordinates(FBtrxxxx instead of FBgnxxxxx).
I am confused now. what else do i need to do to process the matrix of raw counts i just got before importing it into R packages like edgeR for DEG analysis?
jason_ARGONAUTE is offline   Reply With Quote
Old 09-24-2013, 06:56 AM   #5
swaraj
Member
 
Location: Naples, Italy

Join Date: Feb 2012
Posts: 50
Default

In step 2 I coverted transcript BED file into exon file. Here the name of each exon is the transcript name which means transcript name is reduddant. In step 5 this same file is converted back into transcript file using R. All reads in exons with same transcript name are summed up. I would not go into using exon names. Just keep the same name (gene name or trasncript name) for each exon of a feature. This makes it quicker in R to sum the counts of all exons.
swaraj is offline   Reply With Quote
Old 09-25-2013, 08:01 AM   #6
jason_ARGONAUTE
Member
 
Location: china

Join Date: Aug 2013
Posts: 14
Default

ya, i know and i do exactly as ur pipeline. might be i mixed the concept of exon and transcript. i thought they were the same thing.(both named "FBtrxxx")
So, this is not the reason why i got almost twice counts of genes as i expected.
Could you figure out any other reason? Thank you!
jason_ARGONAUTE is offline   Reply With Quote
Old 09-26-2013, 10:44 PM   #7
jason_ARGONAUTE
Member
 
Location: china

Join Date: Aug 2013
Posts: 14
Default

Or,
what do i need to do to get raw counts for each gene?
Thank you!
jason_ARGONAUTE is offline   Reply With Quote
Old 09-27-2013, 03:54 AM   #8
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,476
Default

Quote:
Originally Posted by jason_ARGONAUTE View Post
Or,
what do i need to do to get raw counts for each gene?
Thank you!
If you have a GTF file, just use htseq-count.
dpryan is offline   Reply With Quote
Old 09-27-2013, 06:02 AM   #9
jason_ARGONAUTE
Member
 
Location: china

Join Date: Aug 2013
Posts: 14
Default

I hv tried htseq-count, really hard.
But it cannt be installed on our lab server.
Thank you all the same^^
jason_ARGONAUTE is offline   Reply With Quote
Reply

Tags
bam, bed, counts, coveragebed

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 03:21 AM.


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