The \code{subread-buildindex} (\code{buildindex} function in \Rsubread) program builds an index for reference genome by creating a hash table in which keys are 16bp mers (subreads) extracted from the genome and values are their chromosomal locations.
By default, subreads are extracted from the reference genome at a 2bp interval and a gapped index is built.
During mapping, three sets of subreads will be extracted from each read.
A full index or a gapped index can be built for a reference genome.
In a full index, subreads are extracted from every location in the genome.
In a gapped index, subreads are extracted in every three bases in the genome (ie. there is a 2bp gap between two subreads next to each other).
When a full index is used in read mapping, only one set of subreads are extracted from a read.
However three sets of subreads need to be extracted from a read when a gapped index is used for mapping.
The first set starts from the first base of the read, the second set starts from the second base and the third set starts from the third base.
This makes sure that a mapped read can always have a set of subreads that match those stored in the index.
A full index is larger than a gapped index.
However the full index enables faster mapping speed to be achieved.
When a one-block full index is used for mapping, the maximum mapping speed will be achieved.
Size of one-block full index built for the human reference genome (GRCh38) is 17.8 GB.
Options are available to generate index of any size.
The reference sequences should be in FASTA format.
The \code{subread-buildindex} function divides each reference sequence name (which can be found in the header lines) into multiple substrings by using separators including `\code{|}', ` '(space) and `\code{<tab>}', and it uses the first substring as the name for the reference sequence during its index building.
The first substrings must be distinct for different reference sequences (otherwise the index cannot be built).
Note that the starting `\code{>}' character in the header line is not included in the first substrings.
Sequences of reference genomes can be downloaded from NCBI, UCSC or Ensembl databases.
For instance, the latest mouse reference genome GRCm38.p6 (as of writing) can be downloaded from the NCBI database via \href{ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz}{\color{blue}\url{ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz}}.
Sequences of reference genomes can be downloaded from public databases.
For instance, the primary assembly of human genome GRCh38 or mouse genome GRCm38 can be downloaded from the GENCODE database via the following links:\\
Table 1 describes the arguments used by the \code{subread-buildindex} program.
@@ -384,7 +398,7 @@ Arguments in parenthesis in the first column are used by \code{buildindex}.\newl
\hline
Arguments & Description \\
\hline
chr1.fa, chr2.fa, ... \newline (\code{reference}) & Give names of chromosome files. Note that in{\Rsubread}, only a single FASTA file including all reference sequences should be provided.\\
chr1.fa, chr2.fa, ... \newline (\code{reference}) & Give names of chromosome files. Note for{\Rsubread} only a single FASTA file including all reference sequences should be provided.\\
\hline
-B \newline (\code{indexSplit=FALSE}) & Create one block of index. The built index will not be split into multiple pieces. The more blocks an index has, the slower the mapping speed. This option will override `-M' option when it is also provided.\\
\hline
@@ -392,9 +406,9 @@ chr1.fa, chr2.fa, ... \newline (\code{reference}) & Give names of chromosome fil
\hline
-f $<int>$\newline (\code{TH\_subread}) & Specify the threshold for removing uninformative subreads (highly repetitive 16bp mers). Subreads will be excluded from the index if they occur more than threshold number of times in the reference genome. Default value is 100.\\
\hline
-F \newline (\code{gappedIndex=FALSE}) & Build a full index for the reference genome. 16bp mers (subreads) will be extracted from every position of a reference genome. Size of the full index built for mouse genome is 14GB.\\
-F \newline (\code{gappedIndex=FALSE}) & Build a full index for the reference genome. 16bp mers (subreads) will be extracted from every position of a reference genome.\\
\hline
-M $<int>$\newline (\code{memory}) & Specify the size of computer memory(RAM) in megabytes that will be used for alignment of reads, 8000MB by default. If the size of an index built for a reference genome is greater than the `-M' value, this index will be split into multiple blocks and then saved onto the disk. These blocks will be loaded into computer memory sequentially when performing read alignment. A gapped index generated for mouse genome has a size of 5300MB. Note that when generating a gapped index this function itself uses no more than 12GB of memory.\\
-M $<int>$\newline (\code{memory}) & Specify the size of computer memory(RAM) in megabytes that will be used to store the index during read mapping, 8000MB by default. If the index size is greater than the specified value, the index will be split into multiple blocks. Only one block will be loaded into memory at anytime during the read alignment.\\
\hline
-o $<string>$\newline (\code{basename}) & Specify the base name of the index to be created.\\
\hline
@@ -444,11 +458,11 @@ $^{1,2,3}$ -i $<string> \newline (\code{index}) $ & Specify the base name of the
\hline
$^{1,2}$ -I $<int>$\newline (\code{indels}) & Specify the number of INDEL bases allowed in the mapping. 5 by default. Indels of up to 200bp long can be detected.\\
\hline
$^{1,2}$ -m $<int>$\newline (\code{TH1}) & Specify the consensus threshold, which is the minimal number of consensus subreads required for reporting a hit. The consensus subreads are those subreads which vote for the same location in the reference genome for the read. If pair-end read data are provided, at least one of the two reads from the same pair must satisfy this criteria. 3 by default. For \code{sublong}, this is the consensus threshold for mapping a readlet (1 by default). A readlet is a 100bp sequence extracted from a long read. \\
$^{1,2}$ -m $<int>$\newline (\code{TH1}) & Specify the consensus threshold, which is the minimal number of consensus subreads required for reporting a hit. The consensus subreads are those subreads which vote for the same location in the reference genome for the read. If pair-end read data are provided, at least one of the two reads from the same pair must satisfy this criteria. The default value is 3 for \code{subread-align}, or 1 for \code{subjunc} and \code{sublong}.\\
\hline
$^{1,2}$ -M $<int>$\newline (\code{maxMismatches}) & Specify the maximum number of mis-matched bases allowed in the alignment. 3 by default. Mis-matches found in soft-clipped bases are not counted.\\
\hline
$^{1,2}$ -n $<int>$\newline (\code{nsubreads}) & Specify the number of subreads extracted from each read, 10 by default. For \code{sublong}, this is number of subreads (85 by default) extracted from each readlet. A readlet is a 100bp sequence extracted from a long read.\\
$^{1,2}$ -n $<int>$\newline (\code{nsubreads}) & Specify the number of subreads extracted from each read for mapping. The default value is 10 for \code{subread-align}, or 14 for \code{subjunc}. For \code{sublong}, this is number of subreads (85 by default) extracted from each readlet. A readlet is a 100bp sequence extracted from a long read.\\
\hline
$^{1,2,3}$ -o $<string>$\newline (\code{output\_file}) & Give the name of output file. The default output format is BAM. All reads are included in mapping output, including both mapped and unmapped reads, and they are in the same order as in the input file.\\
\hline
@@ -832,12 +846,9 @@ Such tools use model-based approaches to deconvolve reads overlapping with multi
\subsection{Count multi-mapping reads and multi-overlapping reads}
A multi-mapping read is a read that can be equally best mapped to more than one location in the reference genome.
Due to the mapping ambiguity, it is recommended that multi-mapping reads should be excluded from read counting (default behavior of {\featureCounts} program) to produce as accurate counts as possible.
However we do provide users with other counting options for such reads.
Users can specify the `-M' option (set \code{countMultiMappingReads} to \code{TRUE} in \R) to fully count every alignment reported for a multi-mapping read (each alignment carries 1 count), or specify both `-M' and `--fraction' options (set both \code{countMultiMappingReads} and \code{fraction} to \code{TRUE} in \R) to count each alignment fractionally (each alignment carries $1/x$ count where $x$ is the total number of alignments reported for the read).
Note that for multi-mapping reads the counting is performed at the level of individual alignments (not at read level).
A multi-mapping read is a read that maps to more than one location in the reference genome.
There are multiple options for counting such reads.
Users can specify the `-M' option (set \code{countMultiMappingReads} to \code{TRUE} in \R) to fully count every alignment reported for a multi-mapping read (each alignment carries 1 count), or specify both `-M' and `--fraction' options (set both \code{countMultiMappingReads} and \code{fraction} to \code{TRUE} in \R) to count each alignment fractionally (each alignment carries $1/x$ count where $x$ is the total number of alignments reported for the read), or do not count such reads at all (this is the default behavior in SourceForge {\Subread} package; In \R, you need to set \code{countMultiMappingReads} to \code{FALSE}).
A multi-overlapping read is a read that overlaps more than one meta-feature when counting reads at meta-feature level or overlaps more than one feature when counting reads at feature level.
The decision of whether or not to counting these reads is often determined by the experiment type. We recommend that reads or fragments overlapping more than one gene are not counted for RNA-seq experiments, because any single fragment must originate from only one of the target genes but the identity of the true target gene cannot be confidently determined.
@@ -846,7 +857,7 @@ On the other hand, we recommend that multi-overlapping reads or fragments are co
By default, {\featureCounts} does not count multi-overlapping reads.
Users can specify the `-O' option (set \code{allowMultiOverlap} to \code{TRUE} in \R) to fully count them for each overlapping meta-feature/feature (each overlapping meta-feature/feature receives a count of 1 from a read), or specify both `-O' and `--fraction' options (set both \code{allowMultiOverlap} and \code{fraction} to \code{TRUE} in \R) to assign a fractional count to each overlapping meta-feature/feature (each overlapping meta-feature/feature receives a count of $1/y$ from a read where $y$ is the total number of meta-features/features overlapping with the read).
If a read is both multi-mapping and multi-overlapping, then each overlapping meta-feature/feature will receive a fractional count of $1/(x*y)$ when `-O', `-M', and `--fraction' are all specified.
If a read is both multi-mapping and multi-overlapping, then when `-O', `-M', and `--fraction' are all specified each overlapping meta-feature/feature will receive a fractional count of $1/(x*y)$.
Note that each alignment reported for a multi-mapping read is assessed separately for overlapping with multiple meta-features/features.
When multi-mapping reads are reported with primary and secondary alignments and both `-M' and `--primary' are specified, only primary alignments will be considered in counting and secondary alignments will be ignored.
@@ -886,26 +897,30 @@ When counting reads to meta-features (eg. genes) columns `Chr', `Start', `End' a
Column `Length' always contains one single value which is the total number of non-overlapping bases included in a meta-feature (or a feature), regardless of counting at meta-feature level or feature level.
When counting RNA-seq reads to genes, the `Length' column typically contains the total number of non-overlapping bases in exons belonging to the same gene for each gene.
The counting summary includes the total number of reads that are assigned and also the number of reads that are not assigned due to filtering.
Below lists all the filters supported by {\featureCounts}:
The counting summary includes total number of alignments that were successfully assigned and also number of alignments that failed to be assigned due to various filters.
Note that the counting summary includes the number of alignments, not the number of reads.
Number of alignments will be higher than the number of reads when multi-mapping reads are included since each multi-mapping read contains more than one alignment.
Number and percentage of successfully assigned alignments are also shown in featureCounts screen output.
Filters supported by {\featureCounts} can be found in the list below:
\begin{itemize}
\item Unassigned\_Unmapped: reads are reported as unmapped in SAM/BAM input. Note that if the `--primary' option of featureCounts program is specified, the read marked as a primary alignment will be considered for assigning to features.
\item Unassigned\_MappingQuality: mapping quality scores lower than the specified threshold.
\item Unassigned\_Chimera: two reads from the same pair are mapped to different chromosomes or have incorrect orientation.
\item Unassigned\_FragementLength: length of fragment does not satisfy the criteria.
\item Unassigned\_Duplicate: reads marked as duplicate in the FLAG field in SAM/BAM input.
\item Unassigned\_MultiMapping: reads marked as multi-mapping in SAM/BAM input (the `NH' tag is checked by the program).
\item Unassigned\_Secondary: reads marked as second alignment in the FLAG field in SAM/BAM input.
\item Unassigned\_Nonjunction: reads that do not span exons will not be assigned if the `--countSplitAlignmentsOnly' option is specified.
\item Unassigned\_NoFeatures: not overlapping with any features included in the annotation.
\item Unassigned\_Overlapping\_Length: no features/meta-features were found to have the minimum required overlap length.
\item Unassigned\_Ambiguity: overlapping with two or more features (feature-level summarization) or meta-features (meta-feature-level) summarization.
\item Unassigned\_Unmapped: unmapped reads cannot be assigned.
\item Unassigned\_MappingQuality: alignments with a mapping quality score lower than the threshold.
\item Unassigned\_Chimera: two ends in a paired end alignment are located on different chromosomes or have unexpected orientation.
\item Unassigned\_FragementLength: fragment length inferred from paired end alignment does not meet the length criteria.
\item Unassigned\_Duplicate: alignments marked as duplicate (indicated in the FLAG field).
\item Unassigned\_MultiMapping: alignments reported for multi-mapping reads (indicated by `NH' tag).
\item Unassigned\_Secondary: alignments reported as secondary alignments (indicated in the FLAG field).
\item Unassigned\_Nonjunction: alignments that do not span exons .
\item Unassigned\_NoFeatures: alignments that do not overlap any feature.
\item Unassigned\_Overlapping\_Length: alignments that do not overlap any feature (or meta-feature) with the minimum required overlap length.
\item Unassigned\_Ambiguity: alignments that overlap two or more features (feature-level summarization) or meta-features (meta-feature-level summarization).
\end{itemize}
In the counting summary these filters are listed in the same order as they were applied in counting process (see Section~\ref{sec:read_filtering}).
All categories are exclusive to each other, ie no reads are assigned to more than one category.
If a read can be filtered out by more than one filter, it is always assigned to the first filter it encounters.
All categories are exclusive to each other, ie no alignments are assigned to more than one category.
If an alignment can be filtered out by more than one filter, it is always assigned to the first filter it encounters.
\subsection{Program usage}
@@ -948,7 +963,7 @@ read mapping that produced the provided SAM/BAM files. This optional argument ca
\hline
-L \newline (\code{isLongRead}) & Turn on long-read counting mode. This option should be used when counting long reads such as Nanopore or PacBio reads.\\
\hline
-M \newline (\code{countMultiMappingReads}) & If specified, multi-mapping reads/fragments will be counted. The program uses the `NH' tag to find multi-mapping reads. Alignments reported for a multi-mapping read will be counted separately. Each alignment will have\code{1} count or a fractional count if \code{--fraction} is specified. See section ``Count multi-mapping reads and multi-overlapping reads'' for more details.\\
-M \newline (\code{countMultiMappingReads}) & If specified, multi-mapping reads/fragments will be counted. The program uses the `NH' tag to find multi-mapping reads. Each alignment reported for a multi-mapping read will be counted individually. Each alignment will carry\code{1} count or a fractional count (\code{--fraction}). See section ``Count multi-mapping reads and multi-overlapping reads'' for more details.\\
\hline
-o $<string>$& Give the name of the output file. The output file contains the number of reads assigned to each meta-feature (or each feature if \code{-f} is specified). Note that the {\featureCounts} function in {\Rsubread} does not use this parameter. It returns a \code{list} object including read summarization results and other data. \\
$--$primary \newline (\code{primaryOnly}) & If specified, only primary alignments will be counted. Primary and secondary alignments are identified using bit 0x100 in the Flag field of SAM/BAM files. All primary alignments in a dataset will be counted no matter they are from multi-mapping reads or not (ie. `-M' is ignored).\\
\hline
$--$read2pos $<int>$\newline (\code{read2pos}) &The read is reduced to its 5' most base or 3' most base. Read summarization is then performed based on the single base position to which the read is reduced. By default, no read reduction will be performed.\\
$--$read2pos $<int>$\newline (\code{read2pos}) &Read is reduced to its 5' most base or 3' most base. Read summarization is then performed based on the single base position to which the read is reduced. By default no read reduction is performed. If this option is used with \code{$--$readExtension5} or \code{$--$readExtension3}, a read will be extended first and then reduced to a single base. \\
\hline
$--$readExtension5 $<int>$\newline (\code{readExtension5}) & Reads are extended upstream by $<int>$ bases from their 5' end. 0 by default.\\
\hline
@@ -1205,10 +1220,6 @@ Note that you do not have to run {\repair} before running {\featureCounts}.
The {\repair} program uses a novel approach to quickly find reads from the same pair, rather than performing time-consuming sort of read names.
It takes only about half a minute to re-order a location-sorted BAM file including 30 million read pairs.
\section{coverageCount}
Compute the read coverage for each chromosomal location in the genome.
\section{flattenGTF}
Flatten features (eg. exons) provided in a GTF annotation and output the modified annotation to a SAF format annotation.
SUBREADprintf("\nWarning: failed to find the gene identifier attribute in the 9th column of the provided GTF file.\nThe specified gene identifier attribute is '%s'\nThe attributes included in your GTF annotation are '%s'\n\n",gene_id_column,extra_attrs);
SUBREADprintf("\nERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.\nThe specified gene identifier attribute is '%s'.\nAn example of attributes included in your GTF annotation is '%s'.\nThe program has to terminate.\n\n",gene_id_column,extra_attrs);
}
is_GFF_geneid_warned++;
}
@@ -1127,7 +1129,7 @@ int load_features_annotation(char * file_name, int file_type, char * gene_id_col
SUBREADprintf("\nWarning: failed to find the transcript identifier attribute in the 9th column of the provided GTF file.\nThe specified gene identifier attribute is '%s'\nThe attributes included in your GTF annotation are '%s'\n\n",transcript_id_column,extra_attrs);
SUBREADprintf("\nERROR: failed to find the transcript identifier attribute in the 9th column of the provided GTF file.\nThe specified transcript identifier attribute is '%s'.\nAn example of attributes included in your GTF annotation is '%s'.\nThe program has to terminate\n\n",transcript_id_column,extra_attrs);
}
is_GFF_txid_warned++;
}
@@ -1142,6 +1144,12 @@ int load_features_annotation(char * file_name, int file_type, char * gene_id_col
SUBREADprintf("\nERROR: No feature was loaded from the annotation file. Please check if the annotation format was correctly specified, and also if the feature type was correctly specified if the annotation is in the GTF format.\n\n");
