New upstream version 1.6.5+dfsg

\begin{titlepage}

\begin{center}

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

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

\vspace{1 cm}

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

{\centering\large Rsubread v1.34.6/Subread v1.6.5\\}

\vspace{1 cm}

\centering 11 March 2019\\

\centering 18 July 2019\\

\vspace{5 cm}

\Large Wei Shi and Yang Liao\\

\vspace{1 cm}

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 is achieved.

Size of one-block full index built for the human reference genome (GRCh38) is 17.8 GB.

The \code{subread-buildindex} function needs 15 GB of memory to build this index.

Size of a gapped index built for GRCh38 is less than 9 GB and \code{subread-buildindex} needs 5.7 GB of memory to build it.

Options are available to generate index of any size.

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

\newpage

\section{Memory use and speed}

\label{sec:memoryspeed}

\code{subread-buildindex} (\code{buildindex} function in \Rsubread) needs 15GB of memory to build a full index for human/mouse genome. 

With this index, \code{subread-align} (\code{align} in \Rsubread) require 17.8GB of memory for read mapping.

This enables fastest mapping speed, but it is recommended that the full index should be on a unix server due to relatively large memory use.

Mapping rate is $\sim$14 million reads per minute (10 CPU threads) when full index is used.

A gapped index is recommended for use on a personal computer, which typically has 16GB of memory or less. 

\code{subread-buildindex} (\code{buildindex} function in \Rsubread) only needs 5.7GB of memory to build a gapped index for human/mouse genome. 

\code{subread-align} (\code{align} in \Rsubread) needs 8.2GB of memory for mapping with the gapped index.

It takes \code{subread-buildindex} (\code{buildindex} function in \Rsubread) about 40 minutes to build a full index for human/mouse genome, and building a gapped index takes about 15 minutes.

Memory use for index building and read mapping can be further reduced by building a split index using the \code{-B} and \code{-M} options in \code{subread-buildindex} (\code{indexSplit} and \code{memory} options in \code{buildindex} function in \Rsubread).

\section{Mapping quality scores}

{\Subread} and {\Subjunc} aligners determine the final mapping location of each read by taking into account vote number, number of mis-matched bases, number of matched bases and mapping distance between two reads from the same pair (for paired-end reads only) .

$N_{mm}$ is the number of mismatches present in the final reported alignment for the read.

% \[ MQS = \left\{

% \begin{array}{l l}

% (\sum_{i \in b_m} ( 1 - p_i) - \sum_{i \in b_{mm}} (1 - p_i)) \times 60 / L & \quad \text{if uniquely mapped}\\

% & \quad \text{\scriptsize{[MQS is reset to 0 if less than 0]}}\\

% & \\

% 0 & \quad \text{if $>1$ equally best locations found}\\

% \end{array} \right.\]

% where $L$ is the read length, $p_i$ is the base-calling $p$-value for the $i$th base in the read, $b_m$ is the set of locations of matched bases, and $b_{mm}$ is the set of locations of mismatched bases.

% Base-calling p values can be readily computed from the base quality scores.

% Read bases of high sequencing quality have low base-calling p values.

% Read bases that were found to be insertions are treated as matched bases in the MQS calculation.

% The MQS is a read-length normalized value and it is in the range [0, 60).

\section{Mapping output}

Read mapping results for each library will be saved to a BAM or SAM format file.

\end{Rcode}

\section{Index building}

Please refer to Section~\ref{sec:index}.

Same index is used for the mapping of RNA and DNA sequencing reads. 

\section{Local read alignment}

The \code{Subread} and \code{Subjunc} can both be used to map RNA-seq reads to the reference genome.

Output of {\Subjunc} aligner includes a list of discovered exon-exon junction locations and also the complete alignment results for the reads.

Table 2 describes the arguments used by the {\Subjunc} program.\\

\section{Memory use and speed}

Memory use and running time of \code{subread-buildindex} and \code{subread-align} (\code{buildindex} and \code{align} in \Rsubread) are the same as their memory use and running time in the analysis of DNA sequencing data (see Section~\ref{sec:memoryspeed}).

Compared to \code{subread-align} (\code{align} in \Rsubread), \code{subjunc} uses the same amount of memory when a full index is used and it uses slightly more memory (8.8GB of memory for human/mouse data) when a gapped index is used. 

\code{subjunc} is also slightly slower than \code{subread-align}.

\section{Mapping output}

\item It can count reads at feature (eg. exon) or meta-feature (eg. gene) level.

\item Highly flexible in counting multi-mapping and multi-overlapping reads. Such reads can be excluded, fully counted or fractionally counted.

\item It gives users full control on the summarization of paired-end reads, including allowing them to check if both ends are mapped and/or if the fragment length falls within the specified range.

\item Reduce ambuiguity in assigning read pairs by searching features that overlap with both reads from the pair.

\item Reduce ambiguity in assigning read pairs by searching features that overlap with both reads from the pair.

\item It allows users to specify whether chimeric fragments should be counted.

\item Automatically detect input format (SAM or BAM).

\item Automatically sort paired-end reads. Users can provide either location-sorted or name-sorted bams files to featureCounts. Read sorting is implemented on the fly and it only incurs minimal time cost.

Such tools use model-based approaches to deconvolve reads overlapping with multiple transcripts.

\subsection{Count multi-mapping reads and multi-overlapping reads}

A multi-mapping read is a read that maps to more than one location in the reference genome.

\label{sec:read_filtering}

{\featureCounts} implements a variety of read filters to facilitate flexible read counting, which should satisfy the requirement of most downstream analyses.

The order of these filters being applied is as follows (from highest to lowest):

The order of these filters being applied is as following (from first to last):

unmapped

$>$ read type

$>$ singleton

$>$ mapping quality

$>$ chimeric fragment

$>$ fragment length 

$>$ duplication

$>$ duplicate

$>$ multi-mapping

$>$ secondary alignment

$>$ split reads (or nonsplit reads)

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

The `read type' filter removes those reads that have an unexpected read type and also cannot be counted with confidence.

For example, if there are single end reads included in a paired end read dataset (such data can be produced from a read trimming program for instance) and reads are required to be counted in a strand-specific manner, then all the single end reads will be excluded from counting because their strandness cannot be determined.

However if such reads are to be counted in an unstranded manner then all the single end reads will be considered for counting.

\subsection{Read manipulation}

\begin{itemize}

\item Unassigned\_Unmapped: unmapped reads cannot be assigned.

\item Unassigned\_Read\_Type: reads that have an unexpected read type (eg. being a single end read included in a paired end dataset) and also cannot be counted with confidence (eg. due to stranded counting). Such reads are typically generated from a read trimming program.

\item Unassigned\_Singleton: read pairs that have only one end mapped. 

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

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

An unassigned alignment might fall into more than one category as listed above, however it will only be allocated to one category which is the category corresponding to the first filter that filtered this alignment out. 

\subsection{Program usage}

\section{flattenGTF}

Flatten features (eg. exons) provided in a GTF annotation and output the modified annotation to a SAF format annotation.

If overlapping features are found in the GTF annotation, this function can combine them to form a single large feature encompassing all the original features, or chop them into non-overlapping bins. 

\section{promoterRegions}

#include <assert.h>

#include <stdarg.h>

#include <math.h>

//#include <zutil.h>

#include <unistd.h>

#ifdef MACOS

#include <net/if_dl.h>

#include <netinet/in.h>

#include <arpa/inet.h>

#include <stdlib.h>

#include <string.h>

#include <mach/mach.h>

#include <mach/vm_statistics.h>

#else

#include <netinet/in.h>

#include <net/if.h>

#endif

#include <unistd.h>

#include <sys/types.h>

#include <sys/sysinfo.h>

#endif

#include "gene-algorithms.h"

#include "HelperFunctions.h"

size_t get_sys_mem_info(char * keyword){

	FILE * mfp = fopen("/proc/meminfo","r");

	if(mfp==NULL) return -1;

	char linebuf[1000];

	size_t ret = -1;

	while(1){

		char * rret = fgets(linebuf, 999, mfp);

		if(memcmp( keyword, linebuf, strlen(keyword) ) == 0 && strstr(linebuf," kB")) {

			ret=0;

			int ii ,state=0;

			for(ii=strlen(keyword);; ii++){

				//SUBREADprintf("CH[%d] = %d '%c' at state %d\n", ii, linebuf[ii], linebuf[ii], state);

				if(state == 0 && linebuf[ii]==' ') state = 1;

				if(state == 1 && linebuf[ii]!=' ') state = 2;

				if(state == 2 && linebuf[ii]==' ') state = 9999;

				if(state == 2 && !isdigit(linebuf[ii])){

					SUBREADprintf("WRONG MEMORY INFO '%s'\n", linebuf);

					ret = -1;

					break;

				}

				if(state == 2) ret = ret*10 + ( linebuf[ii] - '0' );

				if(state >= 9999) {

					ret *=1024;

					break;

				}

			}

		}

		if(!rret) break;

	}

	fclose(mfp);

	return ret;

}

char * get_short_fname(char * lname){

	char * ret = lname;

}

#endif

int get_free_total_mem(size_t * total, size_t * free_mem){

#ifdef FREEBSD

    return -1;

#endif

#ifdef MACOS

    mach_msg_type_number_t count = HOST_VM_INFO_COUNT;

    vm_statistics_data_t vmstat;

    int page_size = getpagesize();

    if(KERN_SUCCESS != host_statistics(mach_host_self(), HOST_VM_INFO, (host_info_t)&vmstat, &count))

        return -1;

    //printf("PSIZE=%d\nACT=%u; INACT=%u; FREE=%u\n", page_size, vmstat.active_count, vmstat.inactive_count, vmstat.free_count);

	size_t btlen = sizeof(*total);

    if(sysctl( (int[]) { CTL_HW, HW_MEMSIZE }, 2, total, &btlen, NULL, 0)) return -1;

    *free_mem = (vmstat.free_count + vmstat.inactive_count) * 1llu * page_size;

    return 0;

#else

    struct sysinfo sinf;

    sysinfo(&sinf);

    size_t cached_mem = get_sys_mem_info("Cached:");

    if(cached_mem<0)cached_mem=0;

    *free_mem = cached_mem + sinf.bufferram+sinf.freeram;

    *total = sinf.totalram;

    return 0;

#endif

}

void TNbignum_pow(struct bn* a, struct bn* b, struct bn* c); /* Calculate a^b -- e.g. 2^10 => 1024 */

void TNbignum_isqrt(struct bn* a, struct bn* b);             /* Integer square root -- e.g. isqrt(5) => 2*/

void TNbignum_assign(struct bn* dst, struct bn* src);        /* Copy src into dst -- dst := src */

size_t get_sys_mem_info(char * keyword);

int get_free_total_mem(size_t * total, size_t * free_mem);

#endif

CC = ${CC_EXEC}  ${CCFLAGS}  -fmessage-length=0  -ggdb 

ALL_LIBS= core core-junction core-indel sambam-file sublog gene-algorithms hashtable input-files sorted-hashtable gene-value-index exon-algorithms HelperFunctions interval_merge long-hashtable core-bigtable seek-zlib

ALL_LIBS= core core-junction core-indel sambam-file sublog gene-algorithms hashtable input-files sorted-hashtable gene-value-index exon-algorithms HelperFunctions interval_merge long-hashtable core-bigtable seek-zlib input-blc

ALL_OBJECTS=$(addsuffix .o, ${ALL_LIBS})

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

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