- Count the number of reads in a sam alignment file that map to gff features given a file with aligned sequencing reads and a list of genomic features, a common task is to count how many reads map to each feature. a feature is here an interval (i.e., a range of positions) on a chromosome or a union of such intervals. in the case of rna-seq, the features are typically genes, where each gene is considered here as the union of all its exons. one may also consider each exon as a feature, e.g., in order to check for alternative splicing. for comparative chip-seq, the features might be binding region from a pre-determined list. special care must be taken to decide how to deal with reads that overlap more than one feature. the htseq-count script allows to choose between three modes. of course, if none of these fits your needs, you can write your own script with htseq. see the chapter tour for a step-by-step guide on how to do so. the three overlap resolution modes of htseq-count work as follows. for each position i in the read, a set s(i) is defined as the set of all features overlapping position i. then, consider the set s, which is (with i running through all position within the read) the union of all the sets s(i) for mode union. the intersection of all the sets s(i) for mode intersection-strict. the intersection of all non-empty sets s(i) for mode intersection-nonempty. if s contains precisely one feature, the read is counted for this feature. if it contains more than one feature, the read is counted as ambiguous (and not counted for any features), and if s is empty, the read is counted as no_feature. the following figure illustrates the effect of these three modes: [image]
- Perform simple quality assesment of high-throughput sequencing reads the python script htseq-qa takes a file with sequencing reads (either raw or aligned reads) and produces a pdf file with useful plots to assess the technical quality of a run.