SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BAM files from RNAseq -Alignment (Basespace) not working with DESEq2 in R Rammohan Bioinformatics 6 04-12-2016 11:11 AM
TCGA RNAseq BAM Files fk566938 Bioinformatics 2 12-29-2015 08:27 AM
processing Fastq files from HiSeq2000 single end and RNAseq analysis imsharmanitin RNA Sequencing 1 06-24-2015 09:19 AM
Calling SNPs from RNAseq bam files jmwhitha Bioinformatics 7 10-15-2014 04:16 AM
questions on the three index files for RNAseq analysis wmseq Bioinformatics 6 11-04-2013 08:08 AM

Reply
 
Thread Tools
Old 07-02-2017, 04:42 PM   #1
NDUFB11
Junior Member
 
Location: Italy

Join Date: Jul 2017
Posts: 7
Default RNAseq analysis, and bam files

Hello everyone,
I m handling for the first time bam files from RNAseq and my work is to perform a differential expression analysis. Looking online, I could see there are many tutorials and different ways to do this analysis. I was able to perform example from bioconductor website but I m not able to process my files.

My question is, what is the first step to do with two bam files?

Some tutorial says to use salmon other to use something else,
some others says I have to convert bam to sam files for other bam is enough ...


I would be more than happy if you guys can give me an idea about how to proceed first
thank you
NDUFB11 is offline   Reply With Quote
Old 07-02-2017, 10:33 PM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,466
Default

1. Run featureCounts to produce counts.
2. Load into DESeq2 or edgeR or limma and follow the appropriate vignette

Salmon is a newer method that requires either fastq files or BAM files of transcriptome alignments. It's unlikely that you have those. There is never a reason any more to make a SAM file, that's just a waste of space.
dpryan is offline   Reply With Quote
Old 07-03-2017, 12:19 PM   #3
NDUFB11
Junior Member
 
Location: Italy

Join Date: Jul 2017
Posts: 7
Default

Hi Devon Ryan,
Thank you for your answer, do you know where I can find a good tutorial that explains how to use featureCounts?

or if I can have the comman line

thanks

Last edited by NDUFB11; 07-03-2017 at 06:07 PM.
NDUFB11 is offline   Reply With Quote
Old 07-03-2017, 07:15 PM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

featureCounts manual.
GenoMax is offline   Reply With Quote
Old 07-04-2017, 01:10 PM   #5
NDUFB11
Junior Member
 
Location: Italy

Join Date: Jul 2017
Posts: 7
Default

Thank you genoMax, doesn t seem to be easy :/
NDUFB11 is offline   Reply With Quote
Old 07-04-2017, 04:22 PM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

It is not that difficult if you have some command like skills. If not, can you do R? There is an R package for subread/featurecounts that could be used instead.

All you basically need to do is this (-s 2 is for stranded libraries, if yours is not then use -s 1):
Code:
featureCounts -s 2 -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_file1.bam mapping_results_file2.bam
This will use the annotation file (in GTF format) to identify regions of interest (in this case exon) and then summarize the counts at the gene level (using gene_id) and write the results out to a file (count.txt) as a matrix. Rows would be your gene_id and columns your samples. First 5 columns in the matrix would be annotation. Last two columns will contain the counts for your samples (you have 2 correct, if there are more then include those files in the command line and you will get extra columns in the matrix file). You don't even need to sort the bam files since featureCounts can do that itself.
GenoMax is offline   Reply With Quote
Old 07-04-2017, 06:25 PM   #7
NDUFB11
Junior Member
 
Location: Italy

Join Date: Jul 2017
Posts: 7
Default

Ok I tell you my story,

I have R, containing the package Rsubread::featureCounts() and I have two bam files which I named Cell1.bam and Cell2.bam.

I m in the directory where this files are

> getwd()
[1] "/Users/name/Desktop/bamfolder"


now, if I copy paste this command line I get this
> featureCounts -s 2 -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_file1.bam mapping_results_file2.bam
Error: unexpected numeric constant in "featureCounts -s 2"
>



This is because upstream I miss something, like the word file1.bam and file2.bam need to be Cell1.bam and Cell2.bam, is this correct?

What is missing I guess is the GTF file, should I take it from ensamble?



Thanks a lot

Last edited by NDUFB11; 07-04-2017 at 06:27 PM.
NDUFB11 is offline   Reply With Quote
Old 07-05-2017, 07:35 AM   #8
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

The command I had posted is for using in a shell. You should be able to figure out the R equivalent from subread vignette.

You would need to find a GTF for the genome build that was used to generate your alignment files (you can't mix/match otherwise the results may be nonsense).
GenoMax is offline   Reply With Quote
Old 07-11-2017, 01:31 PM   #9
NDUFB11
Junior Member
 
Location: Italy

Join Date: Jul 2017
Posts: 7
Default

Dear GenoMax,
I think I did another step forward, I run FutureCounts in R.
I used this command line:

> z <- featureCounts(fls, "hg19", nthreads = 30)

but it s not clear to me, how the program counts the reads just writing hg19 as genome reference on the line, without loading any file. Can you tell me more please?

I got this at the end:

========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
Rsubread 1.26.0

//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 2 BAM files ||
|| S wgEncodeRikenCageGm12878CellPapAlnRep1.bam ||
|| S wgEncodeRikenCageHelas3CellPapAlnRep1.bam ||
|| ||
|| Dir for temp files : . ||
|| Threads : 16 ||
|| Level : meta-feature level ||
|| Paired-end : no ||
|| Strand specific : no ||
|| Multimapping reads : not counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
|| ||
|| Load annotation file /usr/local/lib/R/site-library/Rsubread/annot/hg19 ... ||
|| Features : 225074 ||
|| Meta-features : 25702 ||
|| Chromosomes/contigs : 52 ||
|| ||
|| Process BAM file wgEncodeRikenCageGm12878CellPapAlnRep1.bam... ||
|| Single-end reads are included. ||
|| Assign reads to features... ||
|| Total reads : 19677397 ||
|| Successfully assigned reads : 14167723 (72.0%) ||
|| Running time : 0.15 minutes ||
|| ||
|| Process BAM file wgEncodeRikenCageHelas3CellPapAlnRep1.bam... ||
|| Single-end reads are included. ||
|| Assign reads to features... ||
|| Total reads : 24319886 ||
|| Successfully assigned reads : 18864794 (77.6%) ||
|| Running time : 1.08 minutes ||
|| ||
|| Read assignment finished. ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//

Last edited by NDUFB11; 07-11-2017 at 01:34 PM.
NDUFB11 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 06:51 PM.


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