New upstream version 1.6.4+dfsg

\begin{center}

{\Huge\bf Subread/Rsubread Users Guide}\\

\vspace{1 cm}

{\centering\large Subread v1.6.3/Rsubread v1.30.9\\}

{\centering\large Subread v1.6.4/Rsubread v1.32.4\\}

\vspace{1 cm}

\centering 8 October 2018\\

\centering 11 March 2019\\

\vspace{5 cm}

\Large Wei Shi and Yang Liao\\

\vspace{1 cm}

The University of Melbourne\\

Melbourne, Australia\\}

\vspace{7 cm}

\centering Copyright \small{\copyright}  2011 - 2018\\

\centering Copyright \small{\copyright}  2011 - 2019\\

\end{center}

\end{titlepage}

Also included in this software suite is a very efficient SNP caller -- {\ExactSNP}.

{\ExactSNP} measures local background noise for each candidate SNP and then uses that information to accurately call SNPs.

These software programs support a variety of sequencing platforms including Illumina GA/HiSeq, ABI SOLiD, Life Science 454, Helicos Heliscope and Ion Torrent. They are released in two packages -- SourceForge \emph{Subread} package and Bioconductor \emph{Rsubread} package.

These software programs support a variety of sequencing platforms. 

They are released in two packages -- SourceForge \emph{Subread} package and Bioconductor \emph{Rsubread} package\cite{Rsubread}.

\chapter{Preliminaries}

\section{Citation}

If you use {\Subread} or {\Subjunc} aligners, please cite:

If you use {\Rsubread}, you can cite:

\begin{quote}

Liao Y, Smyth GK and Shi W (2013). The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. \emph{Nucleic Acids Research}, 41(10):e108.

Liao Y, Smyth GK and Shi W (2019). The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. \emph{Nucleic Acids Research}, doi: 10.1093/nar/gkz114. [Epub ahead of print].

\\

{\color{blue}{\url{http://www.ncbi.nlm.nih.gov/pubmed/23558742}} }

{\color{blue}{\url{http://www.ncbi.nlm.nih.gov/pubmed/30783653}} }

\end{quote}

{\noindent If you use \featureCounts, please cite:}

If you use \featureCounts, you can cite:

\begin{quote}

Liao Y, Smyth GK and Shi W (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. \emph{Bioinformatics}, 30(7):923-30.

\\

{\color{blue}{\url{http://www.ncbi.nlm.nih.gov/pubmed/24227677}}}

\end{quote}

If you use {\Subread} or {\Subjunc} aligners, you can cite:

\begin{quote}

Liao Y, Smyth GK and Shi W (2013). The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. \emph{Nucleic Acids Research}, 41(10):e108.

\\

{\color{blue}{\url{http://www.ncbi.nlm.nih.gov/pubmed/23558742}} }

\end{quote}

\section{Download and installation}

\subsection{SourceForge {\Subread} package}

\subsection{Install Bioconductor {\Rsubread} package}

{\R} software needs to be installed on my computer before you can install this package.

Launch {\R} and issue the following command to install {\Rsubread}:

\begin{Rcode}

if (!requireNamespace("BiocManager", quietly = TRUE))

    install.packages("BiocManager")

BiocManager::install("Rsubread")

\end{Rcode}

Alternatively, you may download the {\Rsubread} source package directly from {\color{blue}{\url{http://bioconductor.org/packages/release/bioc/html/Rsubread.html}} } and install it from source.\\

To install {\Rsubread} on a Windows computer, the R software has to be launched from Linux Bash Shell (Windows 10) or a unix-like environment such as cygwin.

\subsection{Install SourceForge {\Subread} package}

\subsubsection{Installation from a binary distribution}

This is the easiest way to install the {\Subread} package onto your computer.

Download a {\Subread} binary distribution that suits your oprating system, from the SourceForge website {\color{blue}{\url{http://subread.sourceforge.net}}}. The operating systems currently being supported include multiple variants of Linux (Debian, Ubuntu, Fedora and Cent OS) and Mac OS X. Both 64-bit and 32-bit machines are supported. The executables can be found in the `bin' diretory of the binary package.

This is the easiest way to install the SourceForge {\Subread} package.

Binary distributions can be downloaded for Linux and Mac OS X from {\color{blue}{\url{http://subread.sourceforge.net}}}.

The Linux binary distribution can be run on multiple variants of Linux including Debian, Ubuntu, Fedora and Cent OS.

To install {\Subread} package for other operating systems such as FreeBSD and Solaris, you will have to install them for the source.

For Windows 10, the Linux binary distribution can be directly installed on the Linux Bash Shell. 

If a unix-like environment such as cygwin is installed on a Windows computer, the Linux binary distribution can also be installed in that environment.

To install {\Subread} package on FreeBSD or Solaris, you will have to install them for source.

\subsubsection{Installation from the source package}

\code{make -f Makefile.SunOS}\\

To install it on a Windows computer, you will need to firstly install a unix-like environment such as cygwin and then install the {\Subread} package.\\

If you use a Windows computer, you can install it on the Linux Bash Shell (Windows 10) or install it on a unix-like environment such as cygwin.

The command for installation is the same as that for the installation on a Linux operating system.\\

A new directory called \code{bin} will be created under the home directory of the software package, and the executables generated from the compilation are saved to that directory.

To enable easy access to these executables, you may copy them to a system directory such as \code{/usr/bin} or add the path to them to your search path (your search path is usually specified in the environment variable \code{`PATH'}).

\subsection{Bioconductor {\Rsubread} package}

You have to get {\R} installed on my computer to install this package.

Lauch an {\R} session and issue the following command to install it:

\begin{Rcode}

source("http://bioconductor.org/biocLite.R")

biocLite("Rsubread")

\end{Rcode}

Alternatively, you may download the {\Rsubread} source package directly from {\color{blue}{\url{http://bioconductor.org/packages/release/bioc/html/Rsubread.html}} } and install it to your {\R} from the source.

\section{How to get help}

Bioconductor support site (\url{https://support.bioconductor.org/}) or Google Subread group (\url{https://groups.google.com/forum/#!forum/subread}) are the best place to post questions or make suggestions.

{\noindent\bf Step 1: Build an index}\\

\noindent Build a base-space index (default). You can provide a list of FASTA files or a single FASTA file including all the reference sequences.\\

\noindent Build a base-space index (default). You can provide a list of FASTA files or a single FASTA file including all the reference sequences. The files can be gzipped.\\

\code{subread-buildindex -o my\_index chr1.fa chr2.fa ...}\\

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.

When a one-block full index is used for mapping, the maximum mapping speed is 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.

In \Rsubread, a one-block full index is built by default.

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.

\hline

Arguments & Description \\

\hline

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.\\

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. The files can be gzipped.\\

\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

{\noindent\bf Step 1: Building an index}\\

\noindent The following command can be used to build a base-space index.

You can provide a list of FASTA files or a single FASTA file including all the reference sequences.\\

You can provide a list of FASTA files or a single FASTA file including all the reference sequences.

The files can be gzipped.\\

\code{subread-buildindex -o my\_index chr1.fa chr2.fa ...}\\

{\noindent\bf Step 1: Building an index}\\

\noindent To build the index, you must provide a single FASTA file (eg. ``genome.fa'') which includes all the reference sequences.

The FASTA file can be a gzipped file.

\begin{Rcode}

library(Rsubread)

$>$ duplication

$>$ multi-mapping

$>$ secondary alignment

$>$ split reads

$>$ split reads (or nonsplit reads)

$>$ no overlapping features

$>$ overlapping length

$>$ assignment ambiguity.

Number of reads that were excluded from counting by each filter is reported in the program output, in addition to the reported read counts (see Section~\ref{sec:program_output}). 

\subsection{Read manipulation}

Reads can be shifted (\code{--readShiftType} and \code{--readShiftSize}), extended (\code{--readExtension5} and \code{--readExtension3}) or reduced to an end base (\code{--read2pos}), before being assigned to features/meta-features.

These read manipulations are carried out by {\featureCounts} in the following order: shift $>$ extension $>$ reduction. 

\subsection{Program output}

\label{sec:program_output}

\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\_Split (or Unassigned\_NonSplit): alignments that contain junctions (or do not contain junctions).

\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).

-g $<string>$ \newline (\code{GTF.attrType}) & Specify the attribute type used to group features (eg. exons) into meta-features (eg. genes) when GTF annotation is provided. `gene\_id' by default. This attribute type is usually the gene identifier. This argument is useful for the meta-feature level summarization.\\

\hline

-G $<string>$ \newline (\code{genome}) & Provide the name of a FASTA-format file that contains the reference sequences used in

read mapping that produced the provided SAM/BAM files. This optional argument can be used with '-J' option to improve read counting for junctions.\\

read mapping that produced the provided SAM/BAM files. This optional argument can be used with `-J' option to improve read counting for junctions.\\

\hline

-J \newline (\code{juncCounts}) & Count the number of reads supporting each exon-exon junction. Junctions are identified from those exon-spanning reads (containing `N' in CIGAR string) in input data. The output result includes names of primary and secondary genes that overlap at least one of the two splice sites of a junction. Only one primary gene is reported, but there might be more than one secondary gene reported. Secondary genes do not overlap more splice sites than the primary gene. When the primary and secondary genes overlap same number of splice sites, the gene with the smallest leftmost base position is selected as the primary gene. Also included in the output result are the position information for the left splice site (`Site1') and the right splice site (`Site2') of a junction. These include chromosome name, coordinate and strand of the splice site. In the last columns of the output, number of supporting reads is provided for each junction for each library.\\

-J \newline (\code{juncCounts}) & Count the number of reads supporting each exon-exon junction. Junctions are identified from those exon-spanning reads (containing `N' in CIGAR string) in input data (note that options `--splitOnly' and `--nonSplitOnly' are not considered by this parameter). The output result includes names of primary and secondary genes that overlap at least one of the two splice sites of a junction. Only one primary gene is reported, but there might be more than one secondary gene reported. Secondary genes do not overlap more splice sites than the primary gene. When the primary and secondary genes overlap same number of splice sites, the gene with the smallest leftmost base position is selected as the primary gene. Also included in the output result are the position information for the left splice site (`Site1') and the right splice site (`Site2') of a junction. These include chromosome name, coordinate and strand of the splice site. In the last columns of the output, number of supporting reads is provided for each junction for each library.\\

\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

\hline

$--$byReadGroup \newline (\code{byReadGroup}) & Count reads by read group. Read group information is identified from the header of BAM/SAM input files and the generated count table will include counts for each group in each library.\\

\hline

$--$countSplit \newline AlignmentsOnly \newline (\code{splitOnly}) & If specified, only split alignments (CIGAR strings contain letter `N') will be counted. All the other alignments will be ignored. An example of split alignments is the exon-spanning reads in RNA-seq data. If exon-spanning reads need to be assigned to all their overlapping exons, `-f' and `-O' options should be provided as well.\\

\hline

$--$countNonSplit  \newline AlignmentsOnly \newline (\code{nonSplitOnly}) & If specified, only non-split alignments (CIGAR strings do not contain letter `N') will be counted. All the other alignments will be ignored.\\

\hline

$--$donotsort \newline (\code{autosort}) & If specified, paired end reads will not be re-ordered even if reads from the same pair were found not to be next to each other in the input.\\

\hline

$--$extraAttributes $<string>$ \newline (\code{GTF.attrType.extra}) & Extract extra attribute types from the provided GTF annotation and include them in the counting output. These attribute types will not be used to group features. If more than one attribute type is provided they should be separated by comma (in {\Rsubread} {\featureCounts} its value is a character vector).  \\

\hline

$--$nonOverlapFeature $<int>$ \newline (\code{nonOverlapFeature}) & Maximum number of non-overlapping bases in a feature that is allowed in read assignment. No limit is set by default. \\

\hline

$--$nonSplitOnly \newline (\code{nonSplitOnly}) & If specified, only non-split alignments (CIGAR strings do not contain letter `N') will be counted. All the other alignments will be ignored.\\

\hline

$--$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}) & 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. \\

$--$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. Read reduction is performed after read shifting and read extension if they are also specified. \\

\hline

$--$readExtension3 $<int>$ \newline (\code{readExtension3}) & Reads are extended downstream by $<int>$ bases from their 3' end. 0 by default. Negative value is not allowed. Read extension is performed after read shifting but before read reduction.\\

\hline

$--$readExtension5 $<int>$ \newline (\code{readExtension5}) & Reads are extended upstream by $<int>$ bases from their 5' end. 0 by default.\\

$--$readExtension5 $<int>$ \newline (\code{readExtension5}) & Reads are extended upstream by $<int>$ bases from their 5' end. 0 by default. Negative value is not allowed. \\

\hline

$--$readExtension3 $<int>$ \newline (\code{readExtension3}) & Reads are extended downstream by $<int>$ bases from their 3' end. 0 by default.\\

$--$readShiftSize $<int>$ \newline (\code{readShiftSize}) & Reads are shifted by $<int>$ bases. 0 by default. Negative value is not allowed. \\

\hline

$--$readShiftType $<string>$ \newline (\code{readShiftType}) & Specify the direction in which reads are being shifted. Possible values include \code{upstream}, \code{downstream}, \code{left} and \code{right}.  \code{upstream} by default. Read shifting is performed before read extension or reduction. \\

\hline

$--$Rpath $<string>$ \newline (\code{reportReadsPath}) & Specify a directory to save the detailed assignment results. If unspecified, the directory where counting results are saved is used. See `-R' option for obtaining detailed assignment results for reads.\\

\hline

$--$splitOnly \newline (\code{splitOnly}) & If specified, only split alignments (CIGAR strings contain letter `N') will be counted. All the other alignments will be ignored. An example of split alignments is the exon-spanning reads in RNA-seq data. If exon-spanning reads need to be assigned to all their overlapping exons, `-f' and `-O' options should be provided as well.\\

\hline

$--$tmpDir $<string>$ \newline (\code{tmpDir}) & Directory under which intermediate files are saved (later removed). By default, intermediate files will be saved to the directory specified in `-o' argument (In \R, intermediate files are saved to the current working directory by default).\\

\hline

$--$verbose \newline (\code{verbose}) & Output verbose information for debugging such as unmatched chromosomes/contigs between reads and annotation.\\

\section{removeDup}

Remove duplicated reads from a SAM file.

Remove duplicated reads from a SAM/BAM file. 

In {\Rsubread} this function is called \textsf{removeDupReads}.

\section{subread-fullscan}

ALL_H=$(addsuffix .h, ${ALL_LIBS})

ALL_C=$(addsuffix .c, ${ALL_LIBS})

all: detectionCall sublong repair txUnique featureCounts removeDup exactSNP subread-buildindex subindel  subread-align subjunc qualityScores subread-fullscan propmapped flattenGTF # samMappedBases mergeVCF  testZlib

all: genRandomReads detectionCall sublong repair txUnique featureCounts removeDup exactSNP subread-buildindex subindel  subread-align subjunc qualityScores subread-fullscan propmapped flattenGTF # samMappedBases mergeVCF  testZlib

	mkdir -p ../bin/utilities

	mv longread-one/LRM longread-one/sublong

	mv longread-one/sublong subread-align subjunc featureCounts  subindel exactSNP subread-buildindex ../bin/

	mv detectionCall repair propmapped qualityScores removeDup subread-fullscan txUnique flattenGTF ../bin/utilities

	mv detectionCall genRandomReads repair propmapped qualityScores removeDup subread-fullscan txUnique flattenGTF ../bin/utilities

	@echo

	@echo "###########################################################"

	@echo "#                                                         #"

	rm -f longread-one/*.o

	cd longread-one && $(MAKE)

genRandomReads: gen_rand_reads.c ${ALL_OBJECTS}

	${CC} -o genRandomReads gen_rand_reads.c ${ALL_OBJECTS} ${LDFLAGS}

flattenGTF: flattenAnnotations.c ${ALL_OBJECTS}

	${CC} -o flattenGTF flattenAnnotations.c  ${ALL_OBJECTS} ${LDFLAGS}

ALL_H=$(addsuffix .h, ${ALL_LIBS})

ALL_C=$(addsuffix .c, ${ALL_LIBS})

all: sublong repair featureCounts removeDup exactSNP subread-buildindex subindel  subread-align subjunc qualityScores subread-fullscan propmapped flattenGTF # globalReassembly testZlib

all: genRandomReads sublong repair featureCounts removeDup exactSNP subread-buildindex subindel  subread-align subjunc qualityScores subread-fullscan propmapped flattenGTF # globalReassembly testZlib

	mkdir -p ../bin/utilities

	mv longread-one/LRM longread-one/sublong

	mv longread-one/sublong subread-align subjunc featureCounts  subindel exactSNP subread-buildindex ../bin/

	mv repair subread-fullscan qualityScores removeDup propmapped flattenGTF ../bin/utilities

	mv repair genRandomReads subread-fullscan qualityScores removeDup propmapped flattenGTF ../bin/utilities

	@echo

	@echo "###########################################################"

	@echo "#                                                         #"

	rm -f longread-one/*.o

	cd longread-one && $(MAKE)

genRandomReads: gen_rand_reads.c ${ALL_OBJECTS}

	${CC} -o genRandomReads gen_rand_reads.c ${ALL_OBJECTS} ${LDFLAGS}

flattenGTF: flattenAnnotations.c ${ALL_OBJECTS}

	${CC} -o flattenGTF flattenAnnotations.c  ${ALL_OBJECTS} ${LDFLAGS}