featurecount
Compute the number of reads mapped to genomic features
Syntax
Description
counts the number of reads in the BAM-formatted or SAM-formatted file
T
= featurecount(GTFfile
,Inputfile
)Inputfile
that map onto genomic features as specified in
the GTF-formatted file GTFfile
. GTFfile
specifies the annotation file. Inputfile
specifies the names of
the BAM or SAM files to consider. The output T
is a table where
rows correspond to features and columns correspond to the input files. The elements
of the table consist of the number of reads mapping to each feature for a given
input file.
[___] = featurecount(___,
uses
additional options specified by one or more Name,Value
)Name,Value
pair
arguments.
Examples
Count the number of reads mapped to genomic features
Count reads from a sample SAM file that map to the features included in a GTF file. By default, featurecount
maps the reads to exons, and summarizes the total number of reads at the gene level.
[t,s] = featurecount('Dmel_BDGP5_nohc.gtf','rnaseq_sample1.sam');
Processing GTF file Dmel_BDGP5_nohc.gtf ... Processing SAM file rnaseq_sample1.sam ... Processing reference chr2L ... Processing reference chr2R ... Processing reference chr3L ... Processing reference chr3R ... Processing reference chr4 ... Processing reference chrX ... Done.
Display the first 10 rows of count data.
t(1:10,:)
ans=10×3 table
ID Reference rnaseq_sample1
_______________ _________ ______________
{'FBgn0002121'} {'chr2L'} 9
{'FBgn0067779'} {'chr2L'} 2
{'FBgn0005278'} {'chr2L'} 4
{'FBgn0031220'} {'chr2L'} 4
{'FBgn0025683'} {'chr2L'} 13
{'FBgn0053635'} {'chr2L'} 2
{'FBgn0016977'} {'chr2L'} 22
{'FBgn0086902'} {'chr2L'} 27
{'FBgn0031245'} {'chr2L'} 2
{'FBgn0024352'} {'chr2L'} 2
The ID
column contains the names of features (genes in this example). The Reference
column lists the names of reference sequences for the features. The third column contains the total number of reads mapped to each feature for a given SAM file, that is, rnaseq_sample1.sam. By default, the table shows only those features (rows) and SAM files (columns) with non-zero read counts. Set 'ShowZeroCounts' to true to include those rows and columns with all zero counts in the output table.
s
contains the summary statistics of assigned and unassigned reads from each SAM file. For instance, the TotalEntries
row indicates the total number of alignment records from the given SAM file, and the Assigned
row includes the number of reads that are assigned to features in the GTF file. For details about each row, refer to the Output Arguments section of the reference page.
s
s=9×1 table
rnaseq_sample1
______________
TotalEntries 33354
Assigned 16399
Unassigned_ambiguous 167
Unassigned_filtered 0
Unassigned_lowMappingQuality 0
Unassigned_multiMapped 0
Unassigned_noFeature 16788
Unassigned_supplementary 0
Unassigned_unmapped 0
Count reads without any summarization and disable displaying the progress messages.
[t2,s2] = featurecount('Dmel_BDGP5_nohc.gtf','rnaseq_sample1.sam', ... 'Summarization',false,'Verbose',false);
Notice the ID column of the output table now reports the feature attribute followed by the start and stop positions of each feature, separated by underscores.
t2(1:10,:)
ans=10×3 table
ID Reference rnaseq_sample1
_____________________________ _________ ______________
{'FBgn0002121_12286_12928' } {'chr2L'} 3
{'FBgn0002121_13683_14874' } {'chr2L'} 1
{'FBgn0002121_14933_15711' } {'chr2L'} 3
{'FBgn0067779_67044_67507' } {'chr2L'} 2
{'FBgn0005278_108588_108809'} {'chr2L'} 1
{'FBgn0005278_110755_110877'} {'chr2L'} 1
{'FBgn0005278_112690_113369'} {'chr2L'} 1
{'FBgn0031220_117079_117759'} {'chr2L'} 2
{'FBgn0031220_118361_118874'} {'chr2L'} 1
{'FBgn0031220_118931_119076'} {'chr2L'} 1
You can choose how to assign a read to a particular feature when the read overlaps with multiple features by setting the 'OverlapMethod' option. For instance, if you want to count a read only if it fully overlaps a feature, use the 'full' option.
[tFull, sFull] = featurecount('Dmel_BDGP5_nohc.gtf','rnaseq_sample1.sam', ... 'OverlapMethod','full','Verbose',false);
If you have paired-end data, you can count reads as fragments.
[tFrag,sFrag] = featurecount('Dmel_BDGP5_nohc.gtf','rnaseq_sample1.sam', ... 'CountFragments',true,'Verbose',false);
You can also count fragments from multiple SAM files.
[t2,s2] = featurecount('Dmel_BDGP5_nohc.gtf',... {'rnaseq_sample1.sam','rnaseq_sample2.sam'},'CountFragments',true, ... 'Verbose',false);
Use the following options to count paired-end reads where at least one of the read mates are above a certain mapping quality threshold.
[t3,s3] = featurecount('Dmel_BDGP5_nohc.gtf',... 'rnaseq_sample1.sam','CountFragments',true,'MinMappingQuality',20, ... 'Verbose',false);
If the reads come from any strand-specific assay, you can specify such strand specificity during counting. For instance, if the protocol is stranded, the strand of the feature is compared with the strand of the read. Then only those reads that have the same strand as the overlapped feature are counted.
[t4,s4] = featurecount('Dmel_BDGP5_nohc.gtf',... 'rnaseq_sample1.sam','StrandSpecificity','stranded','Verbose',false);
Input Arguments
GTFfile
— GTF-formatted file name
character vector | string
GTF-formatted file name, specified as a character vector or string.
Example: 'Dmel_BDGP5_nohc.gtf'
Inputfile
— BAM-formatted or SAM-formatted file name
character vector | string | string vector | cell array of character vectors
BAM-formatted or SAM-formatted file name, specified as a character vector, string, string vector, or cell array of character vectors.
Example: 'rnaseq_sample1.sam'
Name-Value Arguments
Specify optional pairs of arguments as
Name1=Value1,...,NameN=ValueN
, where Name
is
the argument name and Value
is the corresponding value.
Name-value arguments must appear after other arguments, but the order of the
pairs does not matter.
Before R2021a, use commas to separate each name and value, and enclose
Name
in quotes.
Example: 'CountFragments',true
specifies to
count reads as pairs of mates.
Feature
— Feature type
'exon'
(default) | character vector | string
Feature type, specified as a character vector or string. This is used to decide what feature to consider from the GTF file. Default is 'exon'
.
Metafeature
— Attribute type
'gene_id'
(default) | character vector | string
Attribute type, specified as a character vector or string. This is used to decide what attribute to consider from the GTF file for grouping features into metafeatures and summarizing the read count.
Summarization
— Boolean variable indicating whether to summarize at the metafeature level
true
(default) | false
Boolean variable indicating whether to summarize at the metafeature level, specified as true
or false
.
Default is true
, meaning the function groups features into metafeatures and reports the read counts for metafeatures.
Alias
— Name of file containing aliases of reference names
character vector | string
Name of file containing aliases of reference names, specified as a character vector or string. The file must be a tab-delimited file where the first column corresponds to the reference names used in the GTF file, and the second column corresponds to the reference names used in the input file(s). The names are case-sensitive. It is necessary to include only the reference names that are different in the GTF file and the input file. The file must contain only one alias term for any reference listed in the input file. By default, the reference names in the GTF file and those in the input files are assumed to be the same.
CountFragments
— Boolean variable indicating whether to count reads as pairs of mates
false
(default) | true
Boolean variable indicating whether to count reads as fragments, specified as true
or false
. Paired-end reads must have the same ID for the field QNAME
in the input file, and the mutual order of mates is inferred by the appropriate bit in the FLAG
field within the input file. Reads that have no valid mate either because the mate is unmapped or filtered out by input criteria are still counted if they satisfy the overlapping criteria.
Default is false
, that is, the reads are counted as single-end reads, and their pairing information is ignored.
StrandSpecificity
— Strand specificity of sequencing protocol
'unstranded'
(default) | 'stranded'
| 'reverse'
Strand specificity of the sequencing protocol, specified as 'unstranded'
(default), 'stranded'
, or 'reverse'
.
If
'unstranded'
, the strand of the reads (or fragments) is ignored.If
'stranded'
, the strand of the reads (or fragments) is considered, and only those having the same strand as the feature they overlap are counted.If
'reverse'
, the opposite direction of the strand of the reads (or fragments) is considered, and only those having the opposite strand as the feature they overlap are counted.
When counting fragments (paired-end reads), the strand of the first mate is considered as the strand of the whole fragment. The mutual order of mates (first or second) is inferred from the appropriate bit in the FLAG
field of the input file.
MinOverlap
— Minimum number of overlapped bases required
1
(default) | positive integer
Minimum number of overlapped bases required to assign a read to a feature, specified as a positive integer. When counting fragments, the sum of the overlaps from each end is used as the minimum number of overlapped bases.
MinMappingQuality
— Minimum mapping quality for a given read
0
(default) | non-negative integer
Minimum mapping quality for a given read to be considered for counting, specified as a non-negative integer. This corresponds to the MAPQ
field in the input file. If counting fragments, at least one of the read mates must satisfy this criterion in order to be considered for counting.
CountMultiOverlap
— Boolean variable indicating whether to count reads overlapping multiple features
false
(default) | true
Boolean variable indicating whether to count reads overlapping multiple features, specified as true
or false
(default).
If true
, a read (or fragment) overlapping multiple features is counted multiple times. During summarization at the metafeature level, a read (or fragment) is counted only once if it overlaps with multiple features belonging to the same metafeature as long as it does not overlap with other metafeaures.
CountMultiMapped
— Counting option for reads having multiple mapping locations in the input file
'primary'
(default) | 'none'
| 'all'
Counting option for reads having multiple mapping locations in the input file, specified as 'primary'
(default), 'none'
, or 'all'
.
If
'primary'
, only the primary alignment of a multi-mapped read is considered. The appropriate bit in the input file is used to identify primary alignments.If
'none'
, all alignments of a multi-mapped read are ignored. The NH tag is used to identify multi-mapped reads.If
'all'
, all alignments of a multi-mapped read are considered and counted multiple times.
BothEndsMapped
— Boolean variable indicating whether a fragment must have both mates mapped
false
(default) | true
Boolean variable indicating whether a fragment must have both mates mapped, specified as true
or false
. Mate mapping information is retrieved from the FLAG
field in the input file. Default is false
.
ProperlyPaired
— Boolean variable indicating whether a fragment must be properly paired
false
(default) | true
Boolean variable indicating whether a fragment must be properly paired, specified as true
or false
. Mate pairing information is retrieved from the FLAG
field in the input file. Default is false
.
ShowZeroCounts
— Boolean variable indicating whether to report features or metafeatures with zero count
false
(default) | true
Boolean variable indicating whether to report features or metafeatures with zero count for every input file in the output table, specified as true
or false
.
Default is false
, that is, only rows with non-zero counts and columns with non-zero counts are included in the output table.
OverlapMethod
— Method to use when assigning a given read to metafeature
'partial'
(default) | 'full'
| 'max'
| 'hits'
Method to use when assigning a given read to metafeature, specified as 'partial'
, 'full'
, 'max'
, or 'hits'
. If 'Summarization'
is set to false
, then the reads are assigned to features, instead of metafeatures, based on the specified method.
In the following table, R refers to a read or fragment, and M refers to a metafeature.
Method | Description |
---|---|
'partial' | R is assigned to M if R overlaps (even partially) only with M. Otherwise R is considered ambiguous. |
'full' | R is assigned to M if R is completely mapped only within M, that is, fully overlapping only M. Otherwise R is considered ambiguous. |
'max' | R is assigned to M if R satisfies the overlapping criteria only with M, or if R satisfies the overlapping criteria with several metafeatures but overlaps fully only with M. |
'hits' | R is assigned to M if R overlaps even partially only M, or if M is the only metafeature with the highest number of features hit by R; otherwise R is considered ambiguous. |
The following schematic diagram and table illustrate the outcome of these methods in conjunction with the 'CountMultiOverlap'
name-value pair argument. In the figure, the read refers to a short-read sequence from an input file, and feature A and feature B refers to features listed in a GTF file.
Each method column lists the feature that the read is assigned to based on the corresponding method. The 'CountMultiOverlap'
column indicates whether this name-value pair is set to true
or false
and if it has any effect in the outcome of each method.
'CountMultiOverlap' | 'partial' | 'full' | 'max' | 'hits' | |
---|---|---|---|---|---|
Case 1 | No effect since the read maps only to one feature (feature A). | feature A | feature A | feature A | feature A |
Case 2 | No effect since the read maps only to one feature (feature A). | feature A | no feature | feature A | feature A |
Case 3 | No effect since the read maps only to one feature (feature A). | feature A | no feature | feature A | feature A |
Case 4 | No effect since the read maps only to one feature (feature A). | feature A | feature A | feature A | feature A |
Case 5 | false | ambiguous | feature A | feature A | ambiguous |
true | feature A, feature B | feature A | feature A | feature A, feature B | |
Case 6 | false | ambiguous | ambiguous | ambiguous | ambiguous |
true | feature A, feature B | feature A, feature B | feature A, feature B | feature A, feature B | |
Case 7 | false | Ambiguous | feature A | feature A | feature A |
true | feature A, feature B | feature A | feature A | feature A |
no feature means that the read is not assigned to any feature. If you have specified the second output table S
, its Unassigned_noFeature
row is incremented by one for such occurrence. ambiguous means that the read is not assigned to any feature since it satisfies the overlapping criteria for multiple features, and the Unassigned_ambiguous
row is incremented by one for such occurrence.
UseParallel
— Boolean variable indicating whether to compute in parallel
false
(default) | true
Boolean variable indicating whether to compute in parallel, specified
as true
or false
.
In order to execute the computation in parallel, you must have Parallel Computing Toolbox™. If a MATLAB® parallel pool does not exist, one is automatically created when the auto-creation option is enabled in your parallel preferences. Otherwise, computation runs in the serial mode.
Default is false
, that is, serial mode.
Verbose
— Boolean variable indicating whether to display the progress of computation
true
(default) | false
Boolean variable indicating whether to display the progress of computation, specified as true
or false
.
Output Arguments
T
— Results containing sequence reads mapped to genomic features
table
Results containing sequence reads mapped to genomic features, returned as a table. The rows correspond to features, and columns correspond to the input files. The elements of the table consist of the number of reads mapped to each feature for a given input file. The table also reports the ID of each feature and the reference sequence for the feature.
When 'Summarization'
is set to false, the
ID column of the table reports the metafeature attribute followed
by the start and stop positions of each feature, separated by underscores.
S
— Summary of assigned and unassigned alignment entries
table
Summary of assigned and unassigned alignment entries, returned as a table. Each column of the table corresponds to each input file provided. The table has the following rows:
Row Names | Description |
---|---|
TotalEntries | Number of records (or alignments) in the input file |
Assigned | Number of reads or fragments that were assigned to features |
Unassigned_ambiguous | Number of unassigned reads or fragments overlapping multiple features or metafeatures |
Unassigned_filtered | Number of alignment records filtered by input criteria |
Unassigned_lowMappingQuality | Number of alignment records filtered out due to low mapping quality |
Unassigned_multiMapped | Number of alignment records not assigned because of corresponding reads mapped to multiple locations |
Unassigned_noFeature | Number of reads or fragments not assigned to any features |
Unassigned_supplementary | Number of alignment records not assigned because they are flagged as supplementary records for chimeric alignments |
Unassigned_unmapped | Number of alignment records not assigned because corresponding reads are unmapped |
Extended Capabilities
Automatic Parallel Support
Accelerate code by automatically running computation in parallel using Parallel Computing Toolbox™.
To run in parallel, set 'UseParallel'
to true
.
For more information, see the 'UseParallel'
name-value pair argument.
Version History
Introduced in R2016a
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)