SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
STAR: ultrafast universal RNA-seq aligner alexdobin Bioinformatics 218 04-02-2018 05:59 PM
Primer design Universal tailed Amplicon paloma84 454 Pyrosequencing 2 11-26-2012 11:33 AM
Program edit an ace file to simultaneously extract read information from all contig? cllorens Bioinformatics 3 06-27-2012 08:20 AM
PubMed: Target-Enrichment Through Amplification of Hairpin-Ligated Universal Targets Newsbot! Literature Watch 0 03-25-2011 06:10 AM
PubMed: Development and quantitative analyses of a universal rRNA-subtraction protoco Newsbot! Literature Watch 1 03-12-2010 05:41 PM

Reply
 
Thread Tools
Old 10-22-2015, 05:42 PM   #141
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,909
Default

Hi Shi,

Brian is on the right track.

Here is a sample of reads

Code:
With v.1.3 tags

XXXXXXX:XXX:XXXXXX:2:1213:18683:18437 1:N:0:        0       chr10   376     3       50M     *       0       0       GTCTGTGTTCTGTGTCAGAGGAGCAGTAGAGCTGAGGCTATCAAGCTATC    BB@BBHEHHCG@<FC@CE@FGH@EGHGFEEHIHFEFE?GEHHHGEHHFHE       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1203:10888:83435 1:N:0:        0       chr10   413     3       50M     *       0       0       CTATCAAGCTATCAGGACATCCCAGATGCCCAAGACTACACTAGTCACCA    DDDBD?HFEE?GHGHEEEFIFIIHGHHHI?HIIIEHIIIIEHHHEHFHHI       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1114:19668:64532 1:N:0:        0       chr10   938     2       5M3I42M *       0       0       CCACTGACACCATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGT    DDDDDIIIIIGIIIIIIIIIHIIIIIIIIHIIIHIIIIIIIIIIIHIIIF       XT:A:R  NM:i:4  AM:i:2
XXXXXXX:XXX:XXXXXX:2:1103:21004:71100 1:N:0:        0       chr10   946     3       50M     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    0<D0@@FHHIIIIIE?F@GHE?H@GHHIIIEGH?HFE@GHHHH?ECHHII       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:2106:11030:47701 1:N:0:        0       chr10   946     3       50M     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    DDDDDHIFIIIHIIIIIIIIHHEHEHIIGHHIIIIHIIHHIIIHHHHHHE       XT:A:R  NM:i:0  AM:i:3

With v.1.4 tags

XXXXXXX:XXX:XXXXXX:2:1213:18683:18437 1:N:0:        0       chr10   376     3       50=     *       0       0       GTCTGTGTTCTGTGTCAGAGGAGCAGTAGAGCTGAGGCTATCAAGCTATC    BB@BBHEHHCG@<FC@CE@FGH@EGHGFEEHIHFEFE?GEHHHGEHHFHE       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1203:10888:83435 1:N:0:        0       chr10   413     3       50=     *       0       0       CTATCAAGCTATCAGGACATCCCAGATGCCCAAGACTACACTAGTCACCA    DDDBD?HFEE?GHGHEEEFIFIIHGHHHI?HIIIEHIIIIEHHHEHFHHI       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:1114:19668:64532 1:N:0:        0       chr10   938     2       2=1X2=3I42=     *       0       0       CCACTGACACCATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTC
AGT     DDDDDIIIIIGIIIIIIIIIHIIIIIIIIHIIIHIIIIIIIIIIIHIIIF      XT:A:R  NM:i:4  AM:i:2
XXXXXXX:XXX:XXXXXX:2:1103:21004:71100 1:N:0:        0       chr10   946     3       50=     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    0<D0@@FHHIIIIIE?F@GHE?H@GHHIIIEGH?HFE@GHHHH?ECHHII       XT:A:R  NM:i:0  AM:i:3
XXXXXXX:XXX:XXXXXX:2:2106:11030:47701 1:N:0:        0       chr10   946     3       50=     *       0       0       ATTCTTACTGTCATCCTCTGAAGCCAGGTCCTGCTCAGTGTTGGCATAGG    DDDDDHIFIIIHIIIIIIIIHHEHEHIIGHHIIIIHIIHHIIIHHHHHHE       XT:A:R  NM:i:0  AM:i:3
Edit: Just saw the last post from Wei. Will look forward to a future release of subread.

Edit2: There was an error in the copy/paste of the data in original post. I have replaced the original with correct copy.

Last edited by GenoMax; 10-23-2015 at 11:03 AM.
GenoMax is offline   Reply With Quote
Old 10-22-2015, 08:19 PM   #142
shi
Wei Shi
 
Location: Australia

Join Date: Feb 2010
Posts: 235
Default

Hi GenoMax,

Thanks for providing the sample reads. However, we noticed that there are missing columns in your mapping results (either PNEXT or TLEN or both). This might cause problems for read counting by featureCounts and lead to incorrect counts. For example, featureCounts might not be able to find NH tags which are used to identify multi-mapping reads.
shi is offline   Reply With Quote
Old 10-23-2015, 10:42 AM   #143
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

That's really strange. It is, indeed, not a valid sam file. I just did some mapping and generated this:

Code:
@HD	VN:1.4	SO:unsorted
@SQ	SN:ecoli_K12	LN:4639675
@PG	ID:BBMap	PN:BBMap	VN:35.51	CL:java -ea -Xmx1g align2.BBMap in=80x_perfect.fa.gz reads=1000 int=f out=xx2.sam ow
400	16	ecoli_K12	778201	45	150=	*	0	0	TCTTCCGGCAACTGATGGACAGGTCAAATTCCCTGCCTGGTCGCCGTATCTGTGATAATAATTAATTGAATAGTAAAGGAATCATTGAAATGCAACTGAACAAAGTGCTGAAAGGGCTGATGATTGCTCTGCCTGTTATGGCAATTGCGG	*	NM:i:0	AM:i:45
401	0	ecoli_K12	940032	45	150=	*	0	0	GTAATACTACTTTCGAGTGAAAATCTACCTATCTCTTTGATTTTCAAATTATTCGATGTATACAAGCCTATATAGCGAACTGCTATAGAAATAATTACACAATACGGTTTGTTACTGGAATCAATCGTGAGCAAGCTTGAGTGAGCCATT	*	NM:i:0	AM:i:45
402	0	ecoli_K12	543578	45	150=	*	0	0	CATTTCATCGCAAATCGCCCGCATGTCGTTTTCTAACTGTTGGGTGAAATCGCGCAGCACGGCAGCGTCGGTATGACGACAATCAATGGTGAACGTGGTTTTACCCGGCACCACATTTACCGTATTCGGGCGCGGCTCTACTTTGCCAAA	*	NM:i:0	AM:i:45
All of the columns are present. I have no idea how you were able to generate output missing some columns... could it be some kind of copy-paste error? Or did you reformat the data in some way?
Brian Bushnell is offline   Reply With Quote
Old 10-23-2015, 10:57 AM   #144
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,909
Default

Quote:
Originally Posted by Brian Bushnell View Post
That's really strange. It is, indeed, not a valid sam file.
All of the columns are present. I have no idea how you were able to generate output missing some columns... could it be some kind of copy-paste error? Or did you reformat the data in some way?
When I copied and pasted the data last night something must have gone wrong. I have recopied the data in post above. Sorry about that.
GenoMax is offline   Reply With Quote
Old 06-02-2016, 04:02 PM   #145
touchsk
Junior Member
 
Location: San Francisco, CA

Join Date: Aug 2015
Posts: 7
Question Question on Output

Hi,
Thanks for this wonderful tool. We use it routinely.
I had a question regarding the output (.featureCounts file). What is the difference between a read having more than one gene assigned with commas in one line (Read1) and a read assigned to more than one gene listed over multiple lines (Read2)?

Read1 Assigned ENSNNNG00000079745,ENSNNNG00000059028 Total=2

Read2 Unassigned_NoFeatures * *
Read2 Unassigned_NoFeatures * *
Read2 Assigned ENSNNNG00000098385 *
Read2 Assigned ENSNNNG00000101430 *
Read2 Assigned ENSNNNG00000100881 *

The output is a result of dropping the --primary option and including the -M and -O options.

Thanks in advance for your help.
Karthik
touchsk is offline   Reply With Quote
Old 06-05-2016, 09:29 PM   #146
shi
Wei Shi
 
Location: Australia

Join Date: Feb 2010
Posts: 235
Default

Read 1 was reported only once in your mapping result and this read was found to overlap with two genes (ENSNNNG00000079745 and ENSNNNG00000059028) by featureCounts.

Read 2 was reported 5 times in your mapping result (5 rows in your bam/sam file). The first two alignments reported for this read do not overlap with any genes and the other three alignments were found to overlap with genes.

Read 1 was counted two times because of multi-overlapping and Read 2 was counted three times due to multi-mapping.
shi is offline   Reply With Quote
Old 06-06-2016, 11:20 AM   #147
touchsk
Junior Member
 
Location: San Francisco, CA

Join Date: Aug 2015
Posts: 7
Default Thanks!

Thanks for the quick response and clarification provided.
touchsk is offline   Reply With Quote
Old 10-30-2016, 09:24 AM   #148
htetre
Member
 
Location: US

Join Date: Jul 2013
Posts: 28
Default

Hello,

I am trying to use featureCount with a bam file that aligned using Hisat2 with an indexed genome and spliced sites file created from a gtf file. Now I am using the same gtf file in featureCount but am getting 0% Successfully assigned reads.

Here is the code I am using
featureCounts -a gene.gtf -F GTF -b -o countsSb10.txt -t gene -R -g gene_id Sb10.aln-sort.bam

I am not sure what I am doing wrong. These are reads that aligned to the reference genome, my unmapped reads are removed.

Thank you
htetre
htetre is offline   Reply With Quote
Old 10-31-2016, 02:37 PM   #149
shi
Wei Shi
 
Location: Australia

Join Date: Feb 2010
Posts: 235
Default

Are you using the latest version? '-b' is an invalid argument. Could you also provide the assignment summary that can be found in the 'countsSb10.txt.summary' file?
shi is offline   Reply With Quote
Old 05-02-2017, 09:05 AM   #150
TML
Junior Member
 
Location: PNW

Join Date: May 2014
Posts: 1
Default output tagged bamfile

Hi Shi,

Would it be possible to include an option for outputting a tagged bam file (similar to the htseq --samout option) in the next release. Basically a bamfile with a XF:Z: tag that has the name of the assignment there. We're using this currently to allow umi-tools pcr deduplication of RNAseq reads on a per gene basis, and use the FC -R files to add the tag to the bamfiles ourselves but it would a lot faster if featurecounts could just output the tagged bamfiles.

Thanks, Marinus
TML is offline   Reply With Quote
Old 03-07-2019, 01:53 PM   #151
rookie_genomics
Junior Member
 
Location: USA

Join Date: Mar 2019
Posts: 4
Default

Hi,

I am new to RNA sequencing analysis and just finished assigning my reads using featureCounts. I am trying to save the output of featureCounts into a txt file and am having trouble.

The command I used for fc is as follows

counts <- featureCounts("my.bam",annot.ext ="my.gtf",isGTFAnnotationFile=T,GTF.featureType="exon",GTF.attrType="gene_id",nthreads=4,isPairedEnd=T,countMultiMappingReads=T)

And I used this command to tabulate the results

write.table(cbind(counts$annotation[,2:4], counts$counts),"sample_featureCounts.txt",quote=F,sep=" ",row.names=F)

My output file has only 4 columns

Chromosome(Ensembl ID) Start End Counts

I would like to add gene name to the output to identify counts better.

How can I do that? What command should I use?
rookie_genomics is offline   Reply With Quote
Old 03-07-2019, 02:39 PM   #152
yangliao
Member
 
Location: Melboune, VIC

Join Date: May 2013
Posts: 13
Default

Quote:
Originally Posted by rookie_genomics View Post
Code:
write.table(cbind(counts$annotation[,2:4], counts$counts),"sample_featureCounts.txt",quote=F,sep=" ",row.names=F)
You saved the 2nd to the 4th columns in the annotation data frame into the file. These three columns are the chromosome names, start and end locations, as you had in the output.

The gene identity is in the first column of annotation, so you can use
Code:
write.table(cbind(counts$annotation[,1:4], counts$counts),"sample_featureCounts.txt",quote=F,sep=" ",row.names=F)
yangliao is offline   Reply With Quote
Old 04-09-2019, 06:34 AM   #153
Guillaume
Junior Member
 
Location: Marseille, France

Join Date: Mar 2010
Posts: 2
Default COunt inconsistency

Hello

I have problems getting correct counts for a stranded RNA-seq study using hisat2 version 2.1.0 and featureCounts v1.6.4, and I would appreciate any help.

You have below an image of the results of the mapping by Hisat2 (--rna-strandness RF) of a small test subset of my data

For example, the second last gene (AFBG1_15566.1) is only covered by transcripts of the same orientation (blue reads).

However, when I process the same BAM file with featureCounts, it finds
with the flag -s 1:

AFBG1_15566 scf7180000002748;scf7180000002748 13827;16161 16084;16780 +;+ 2878 2383

with the flag -s 2:
AFBG1_15566 scf7180000002748;scf7180000002748 13827;16161 16084;16780 +;+ 2878 2257


I can't understand why featureCounts finds about the same counts (2383 and 2257) in both orientations.


if that can be useful to find the solution, I have put the script and data that allow to get these results at this address: https://nuage.osupytheas.fr/s/D5S4stT9aDLYEsD

Thank you !


Last edited by Guillaume; 04-09-2019 at 06:37 AM.
Guillaume is offline   Reply With Quote
Old 04-09-2019, 02:12 PM   #154
yangliao
Member
 
Location: Melboune, VIC

Join Date: May 2013
Posts: 13
Default

Hi Guillaume,

I found that you did not use the "-p" option in your script to run featureCounts. This means that each single read, not read-pair, was assigned to the genes. FeatureCounts only flips the strand of the second read when the "-p" option is specified; otherwise it simply looks for the 0x10 FLAG in the alignment for matching the strands of the gene and the alignment. Half of your single reads (R2s) were from the positive strand, while the other half (R1s) were from negative strand, hence the AFBG1_15566 gene always has counts no matter you used "-s 1" or "-s 2".

I changed your script by adding the "-p" option to featureCounts, flipping all your R2s to the negative strand. Now using "-s 1" has zero count for AFBG1_15566, but using "-s 2" has all counts for AFBG1_15566.

Cheers,
Yang
yangliao is offline   Reply With Quote
Old 04-10-2019, 04:38 AM   #155
Guillaume
Junior Member
 
Location: Marseille, France

Join Date: Mar 2010
Posts: 2
Default

Thank you so much Yang
I should have read the manual more carefully...
Guillaume 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:52 AM.


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