| Safe Haskell | None |
|---|---|
| Language | Haskell2010 |
Bio.Adna
- data DmgStats a = DmgStats {
- basecompo5 :: CompositionStats
- basecompo3 :: CompositionStats
- substs5 :: SubstitutionStats
- substs3 :: SubstitutionStats
- substs5d5 :: SubstitutionStats
- substs3d5 :: SubstitutionStats
- substs5d3 :: SubstitutionStats
- substs3d3 :: SubstitutionStats
- substs5dd :: SubstitutionStats
- substs3dd :: SubstitutionStats
- substs5cpg :: SubstitutionStats
- substs3cpg :: SubstitutionStats
- stats_more :: a
- type CompositionStats = [(Maybe Nucleotide, Vector Int)]
- type SubstitutionStats = [(Subst, Vector Int)]
- addFragType :: BamMeta -> Enumeratee [BamRaw] [(BamRaw, FragType)] m b
- damagePatternsIter :: MonadIO m => Int -> Int -> Iteratee [Alignment] m b -> Iteratee [(BamRec, FragType, Vector Word8, Vector NPair)] m (DmgStats b)
- damagePatternsIterMD :: MonadIO m => Int -> Iteratee [Alignment] m b -> Iteratee [(BamRaw, FragType)] m (DmgStats b)
- damagePatternsIter2Bit :: MonadIO m => Refs -> TwoBitFile -> Int -> Int -> Iteratee [Alignment] m b -> Iteratee [(BamRaw, FragType)] m (DmgStats b)
- alnFromMd :: Vector_Nucs_half Nucleotides -> Vector Cigar -> [MdOp] -> Vector NPair
- data DamageParameters float = DP {
- ssd_sigma :: !float
- ssd_delta :: !float
- ssd_lambda :: !float
- ssd_kappa :: !float
- dsd_sigma :: !float
- dsd_delta :: !float
- dsd_lambda :: !float
- data NewDamageParameters vec float = NDP {}
- data GenDamageParameters vec float
- = UnknownDamage
- | OldDamage (DamageParameters float)
- | NewDamage (NewDamageParameters vec float)
- type DamageModel = Bool -> Int -> Int -> Mat44D
- bang :: Mat44D -> Subst -> Double
- nudge :: MMat44D -> Subst -> Double -> IO ()
- data Alignment = ALN {
- a_sequence :: !(Vector NPair)
- a_fragment_type :: !FragType
- data FragType
- data Subst = Nucleotide :-> Nucleotide
- type NPair = (Nucleotides, Nucleotides)
- noDamage :: DamageModel
- univDamage :: DamageParameters Double -> DamageModel
- empDamage :: NewDamageParameters Vector Double -> DamageModel
- newtype Mat44D = Mat44D (Vector Double)
- newtype MMat44D = MMat44D (IOVector Double)
- scalarMat :: Double -> Mat44D
- complMat :: Mat44D -> Mat44D
- freezeMats :: MMat44D -> MMat44D -> IO Mat44D
- bwa_cal_maxdiff :: Double -> Int -> Int
Documentation
Collected "traditional" statistics:
- Base composition near 5' end and near 3' end. Each consists of
five vectors of counts of A,C,G,T, and everything else.
basecompo5begins withcontextbases to the left of the 5' end,basecompo3ends withcontextbases to the right of the 3' end. - Substitutions. Counted from the reconstructed alignment, once around the 5' end and once around the 3' end. For a total of 2*4*4 different substitutions. Positions where the query has a gap are skipped.
- Substitutions at CpG motifs. Also counted from the reconstructed alignment, and a CpG site is simply the sequence CG in the reference. Gaps may confuse that definition, so that CpHpG still counts as CpG, because the H is gapped. That might actually be desirable.
- Conditional substitutions. The 5' and 3' ends count as damaged if
the very last position has a C-to-T substitution. With that in
mind,
substs5d5,substs5d3,substs5ddare likesubsts5, but counting only reads where the 5' end is damaged, where the 3' end is damaged, and where both ends are damaged, respectively.
XXX This got kind of ugly. We'll see where this goes...
Constructors
type CompositionStats = [(Maybe Nucleotide, Vector Int)] Source #
addFragType :: BamMeta -> Enumeratee [BamRaw] [(BamRaw, FragType)] m b Source #
damagePatternsIter :: MonadIO m => Int -> Int -> Iteratee [Alignment] m b -> Iteratee [(BamRec, FragType, Vector Word8, Vector NPair)] m (DmgStats b) Source #
Common logic for statistics. The function get_ref_and_aln
reconstructs reference sequence and alignment from a Bam record. It
is expected to construct the alignment with respect to the forwards
strand of the reference; we reverse-complement it if necessary.
damagePatternsIterMD :: MonadIO m => Int -> Iteratee [Alignment] m b -> Iteratee [(BamRaw, FragType)] m (DmgStats b) Source #
Enumeratee (almost) that computes some statistics from plain BAM
with a valid MD field. The Alignment is also reconstructed and
passed downstream. The result of any downstream processing is
available in the stats_more field of the result.
- Reconstruct the alignment from CIGAR, SEQ, and MD.
- Filter the alignment to get the reference sequence, accumulate it.
- Accumulate everything over the alignment.
The argument is the amount of context inside desired.
For Complete fragments, we cut the read in the middle, so the 5'
and 3' plots stay clean from each other's influence. Leading and
Trailing fragments count completely towards the appropriate end.
damagePatternsIter2Bit :: MonadIO m => Refs -> TwoBitFile -> Int -> Int -> Iteratee [Alignment] m b -> Iteratee [(BamRaw, FragType)] m (DmgStats b) Source #
Enumeratee (almost) that computes some statistics from plain BAM
(no MD field needed) and a 2bit file. The Alignment is also
reconstructed and passed downstream. The result of any downstream
processing is available in the stats_more field of the result.
- Get the reference sequence including both contexts once. If this includes invalid sequence (negative coordinate), pad suitably.
- Accumulate counts for the valid parts around 5' and 3' ends as appropriate from flags and config.
- Combine the part that was aligned to (so no context) with the read to reconstruct the alignment.
Arguments are the table of reference names, the 2bit file with the reference, the amount of context outside the alignment desired, and the amount of context inside desired.
For Complete fragments, we cut the read in the middle, so the 5'
and 3' plots stay clean from each other's influence. Leading and
Trailing fragments count completely towards the appropriate end.
alnFromMd :: Vector_Nucs_half Nucleotides -> Vector Cigar -> [MdOp] -> Vector NPair Source #
Reconstructs the alignment from query, cigar, and md. Only positions where the query is not gapped are produced.
data DamageParameters float Source #
Parameters for the universal damage model.
We assume the correct model is either no damage, or single strand damage, or double strand damage. Each of them comes with a probability. It turns out that blending them into one is simply accomplished by multiplying these probabilities onto the deamination probabilities.
For single stranded library prep, only one kind of damage occurs (C
frequency (ssd_sigma) in single stranded parts, and the overhang
length is distributed exponentially with parameter ssd_lambda at
the 5' end and ssd_kappa at the 3' end. (Without UDG treatment,
those will be equal. With UDG, those are much smaller and in fact
don't literally represent overhangs.)
For double stranded library prep, we get C->T damage at the 5' end
and G->A at the 3' end with rate dsd_sigma and both in the interior
with rate dsd_delta. Everything is symmetric, and therefore the
orientation of the aligned read doesn't matter either. Both
overhangs follow a distribution with parameter dsd_lambda.
Constructors
| DP | |
Fields
| |
Instances
| Read float => Read (DamageParameters float) Source # | |
| Show float => Show (DamageParameters float) Source # | |
| Generic (DamageParameters float) Source # | |
| type Rep (DamageParameters float) Source # | |
data NewDamageParameters vec float Source #
Constructors
| NDP | |
Instances
| (Read (vec float), Read float) => Read (NewDamageParameters vec float) Source # | |
| (Show (vec float), Show float) => Show (NewDamageParameters vec float) Source # | |
| Generic (NewDamageParameters vec float) Source # | |
| type Rep (NewDamageParameters vec float) Source # | |
data GenDamageParameters vec float Source #
Constructors
| UnknownDamage | |
| OldDamage (DamageParameters float) | |
| NewDamage (NewDamageParameters vec float) |
Instances
| (Read (vec float), Read float) => Read (GenDamageParameters vec float) Source # | |
| (Show (vec float), Show float) => Show (GenDamageParameters vec float) Source # | |
| Generic (GenDamageParameters vec float) Source # | |
| type Rep (GenDamageParameters vec float) Source # | |
type DamageModel = Bool -> Int -> Int -> Mat44D Source #
A DamageModel is a function that gives substitution matrices for
each position in a read. The DamageModel can depend on whether the
alignment is reversed, the length of the read and the position. (In
practice, we should probably memoize precomputed damage models
somehow.)
bang :: Mat44D -> Subst -> Double infix 8 Source #
Convenience function to access a substitution matrix that has a mnemonic reading.
Constructors
| ALN | |
Fields
| |
Constructors
| Nucleotide :-> Nucleotide infix 9 |
type NPair = (Nucleotides, Nucleotides) Source #
noDamage :: DamageModel Source #
DamageModel for undamaged DNA. The likelihoods follow directly
from the quality score. This needs elaboration to see what to do
with amibiguity codes (even though those haven't actually been
observed in the wild).
Things specific to ancient DNA, e.g. damage models.
For aDNA, we need a substitution probability. We have three options: use an empirically determined PSSM, use an arithmetically defined PSSM based on the Johnson model, use a context sensitive PSSM based on the Johnson model and an alignment. Using Dindel, actual substitutions relative to a called haplotype would be taken into account. Since we're not going to do that, taking alignments into account is difficult, somewhat approximate, and therefore not worth the hassle.
We represent substitution matrices by the type Mat44D. Internally,
this is a vector of packed vectors. Conveniently, each of the packed
vectors represents all transition into the given nucleotide.
freezeMats :: MMat44D -> MMat44D -> IO Mat44D Source #
Adds the two matrices of a mutable substitution model (one for each strand) appropriately, normalizes the result (to make probabilities from pseudo-counts), and freezes that into one immutable matrix. We add a single count everywhere to avoid getting NaNs from bizarre data.
bwa_cal_maxdiff :: Double -> Int -> Int Source #
Number of mismatches allowed by BWA.
bwa_cal_maxdiff thresh len returns the number of mismatches
bwa aln -n $tresh would allow in a read of length len. For
reference, here is the code from BWA that computes it (we assume err
= 0.02, just like BWA):
int bwa_cal_maxdiff(int l, double err, double thres)
{
double elambda = exp(-l * err);
double sum, y = 1.0;
int k, x = 1;
for (k = 1, sum = elambda; k < 1000; ++k) {
y *= l * err;
x *= k;
sum += elambda * y / x;
if (1.0 - sum < thres) return k;
}
return 2;
}
double sum, y = 1.0; int k, x = 1; for (k = 1, sum = elambda; k < 1000; ++k) { y *= l * err; x *= k; sum += elambda * y / x; if (1.0 - sum < thres) return k; } return 2; } @