Question: How to generate a count matrix with featurecounts
1
3
Entering edit mode
6.0 years ago
kamel ▴ 70

Hi, I used featureocunts for quantization but its output is complicated. PLEASE I need to generate a count matrix with featurecounts

     

RNA-Seq rna-seq • 15k views
ADD COMMENT
0
Entering edit mode

If you used featureCounts and got a complicated output, maybe it is already the count matrix? How did you run it? Whyyou think the output you got is not a counts matrix? Can you show a snippet of the output?

ADD REPLY
0
Entering edit mode

maybe I did not explain it well. I want to keep only the gene-id and the number of reads assigned

Here are the results I found:

Geneid  Chr     Start   End     Strand  Length ./bam
ASMNOP1_01G00001       1;1;1;1;1;1;1;1;1;1;1   6089;6365;6728;6865;7041;7417;7727;8116;8595;8658;10466 6096;6676;6812;6972;7358;7676;8063;8523;8609;10408;11000        -;-;-;-;-;-;-;-;-;-;-   4137    9
ASMNOP1_01G00003       1;1;1;1 11394;11702;11793;12130 11649;11736;12083;12168 -;-;-;- 621     0
ASMNOP1_01G00005       1       19069   19668   -       600     0
ASMNOP1_01G00007       1;1;1;1;1       21463;21835;21912;22037;22138   21516;21874;21959;22060;22193   +;+;+;+;+       222     0
ASMNOP1_01G00009       1;1;1;1;1;1;1;1 23519;23690;23942;24147;24382;24550;24][1]
ADD REPLY
0
Entering edit mode

The output of featureCounts is a tab-separated file, you can use cut to get the columns you want. By the way, you know if you input several bams to featureCounts, its output will include one column of counts for each input bam?

ADD REPLY
0
Entering edit mode

I do not know how I can do cat in this situation. would it be possible to give me the command that I can use to have that number of reads of each gene-id.

  I used two bam files and my results fit like this:

Geneid Chr     Start   End     Strand  Length   .library2/bam  ./library1/bam
    ASMNOP1_01G00001 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1 6089; 6365; 6728; 6865; 7041; 7417; 7727; 8116; 8595; 8658; 10466 6096; 6676; 6812; 6972 ; 7358; 7676; 8063; 8523; 8609; 10408; 11000 -; -; -; -; -; -; -; -; -; 4137 9 59
    ASMNOP1_01G00003 1; 1; 1; 1 11394; 11702; 11793; 12130 11649; 11736; 12083; 12168 -; -; -; - 621 0 0
    ASMNOP1_01G00005 1 19069 19668 - 600 0 0\]\[1\]\]\[1\]][1]
ADD REPLY
0
Entering edit mode

Numbers under the columns with .library2/bam and ./library1/bam should be the counts for genes in column 1. (How did you end up with strange sample file names like that?)

ADD REPLY
0
Entering edit mode

here is the command I used: Program:featureCounts v1.6.1; Command: featureCounts" "-T" "12" "-s" "2" "-t" "exon" "-g" "gene_id" "-a" "/annotation.gtf" "-o" "counts.txt" "/library1.bam" "library2.bam"enter code here

ADD REPLY
1
Entering edit mode

While slightly mangled by the quotes you included around each option the command looks ok. Like I said above last two columns in your file have the counts. You can read this matrix file into DESeq2 etc and manipulate it there (or use cut -f1,7-8 your_file > new_file to extract gene_ID and count columns in a new file).

ADD REPLY
0
Entering edit mode

Thank you very much, you helped me. yes I will use DESeq2 so easily.

a quick question for multimapper reads, I had the idea that featureCounts does not count reads but I found that with the option "-M" it's what it allows to count all reads multimappers.

ADD REPLY
3
Entering edit mode

While that is correct you should not count multi-mapping alignments for RNAseq (featureCounts for paired bam files and dealing with multimapping reads? ).

ADD REPLY
1
Entering edit mode

Or if you want to count multi-mappers, do it with a program that distributes the counts properly, such as RSEM or Salmon - however, you need to map to the transcriptome, then.

ADD REPLY
0
Entering edit mode

I can not count the multimappers reads by a mapping on the genome?

ADD REPLY
0
Entering edit mode

Thank you for your reply genomax. I think I need to count multimapper reads because I'm working on a hexaploid genome

ADD REPLY
0
Entering edit mode

Those recommendations are for any genome in general.

ADD REPLY
0
Entering edit mode

cut is not cat.

ADD REPLY
0
Entering edit mode

SorryI wanted to write cut. can you give me the command that I can use

ADD REPLY
1
Entering edit mode

cut -f and the column numbers you want to select, separated by commas (,).

man cut for details.

ADD REPLY
0
Entering edit mode

Thank you very much.

ADD REPLY
4
Entering edit mode
ADD COMMENT

Login before adding your answer.

Traffic: 3013 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6