full
border
#666666
http://www.novocraft.com/wp-content/themes/smartbox-installable/
http://www.novocraft.com/
#0397c9
style1

NovoMethyl – Analyzing Methylation Status

Based on Novomethyl (Beta.9.0)

Introduction

Novomethyl can analyse a set of alignments to identify methylated cytosine’s.

Novomethyl analyses Cytosines using two models

Methylation State In this model the sample is assumed to be homogeneous and we call cytosines as either methylated or unmethylated.  Each base in the genome being sequenced can be one of A, Cm, Cu, Gm, Gu or T and calling methylation status is very similar to calling a consensus sequence or SNPs with samtools pileup except that we have 6 rather than 4 states for each position. Novomethyl adopts the base calling models from progams such as samtools pileup to call methylation status.In this case we will use the basic statistical models from SOAP SNP detection for massively parallel whole-genome resequencing” (2009) Genome Res. published online May 6, 2009 , doi:10.1101/gr.088013.108
Methylation Percentage This model is designed for heterogeneous samples and uses a traditional 4 base model to call the consensus sequence and then for each cytosine it reports the % of bases that were methylated (i.e. Not converted to thyamine).  This model should be used for heterogenous samples such as tumor cells.

In both modes the reference sequence is used to set prior probabilities for the base calls and then the actual base calls are used to calculate posterior probabilities. The models can handle SNPs and heterozygous locii. It’s fully possible that cytosines in the reference are not in present in the sample or vice versa.

 

CT & GA alignments

DNA fragments can come off the Watson or Crick strand of the sample however our reference sequence has only one strand.  Alignments of fragments from the Watson strand align in CT mode (i.e. T’s in the read can align to either a C or T in the reference sequence) and fragments from the Crick strand align in GA mode because they need to be reverse complemented to align to the reference and now Crick strand cytosine’s are the G’s in the reference.

On the Watson strand we represent methylated cytosine as Cm and unmethylated as Cu and Cytosines from the Crick strand are represented as Gm and Gu depending on methylation status.

 

Report Formats

There are two basic report formats, BED format and a native report format that borrows from samtools pileup format with addition of extra fields for methylation status.

The BED format report comes in two flavours:

  1. Track name is Cytosine Context & Methylation State and score is Phred scaled quality of the state
  2. Track name is Cytosine Context and the the score is the Methylation Percentage

 

Method

Four basic steps

  1. Align using Novoalign V2.07.06 or later
  2. Split alignments into separate CT & GA BAM files
  3. Use samtools mpileup to produce a raw pileup from the CT & GA files
  4. Use novomethyl to call methylation status

 

Aligning

Nothing special here, just make sure you use Novoalign V2.07.06 or later as it tags alignments based on the Watson/Crick strand. For best performance use -b2 option. we need to use a recent copy of Novoalign as it adds a tag to the sam records indication if the alignment is in CT or GA mode.

 

Splitting the BAM files

This is necessary as we need to separate alignments into Watson & Crick strands. This is easily done using grep on the sam report from Novoalign, we just grep for lines with tag ZB:Z:CT or ZB:Z:GA.

Example code that processes a set of single end read files:

 

for f in {list of read files} 
do 
    novoalign -d referencegenome.nbx -t 20,3 -f ${f} -F STDFQ -b2 -a -o SAM '@RG\tID:Pooh' > ${f}.sam 2>${f}.log 
    grep -e "^@\|ZB:Z:CT" ${f}.sam | samtools view -uS - | samtools sort - ${f}.CT& 
    grep -e "^@\|ZB:Z:GA" ${f}.sam | samtools view -uS - | samtools sort - ${f}.GA& 
done 
wait

The option -t 20,3 sets the alignment threshold. For Bi-Seq projects we usually set this lower than the default.

The option -a turns on adapter trimming for short fragments. This is especially useful if you’ve gone for a short insert length to try and avoid fragments breaking during bisulphite treatment. For paired end  data sets we’ve tested the default adapters were correct. However we suggest you check what adapters were actually used.

After samtools sort you may optionally remove duplicate fragments with samtools rmdup.

If you had multiple read files the next step is to merge the BAM files from above (but keeping CT & GA files separate)

Skip this step if you only had one file of reads.

samtools merge ct.bam *.CT.bam&
samtools merge ga.bam *.GA.bam&
wait

Samtools mpileup

We now have two BAM files and we can use these with samtools mpileup to produce a raw mpileup file. Note that we are not using mpileup to call consensus or bases, just to pileup the bases.

samtools mpileup -BC 0 -q 30 -f referencegenome.fa ct.bam ga.bam >biseq.mpileup.txt

The mpileup options -BC 0 are required to turn off base quality calibration. If not used, then mpileup will lower the quality of converted cytosine calls.

We’ve also shown the use of a quality limit of 30 for alignments, this is more an intuition than a recommendation.

A couple of important points:

  1. Don’t use an SM tag on the @RG record
  2. The ct.bam file should be first

The mpileup file will have separate columns for Watson & Crick strands (CT & GA mode alignments) and  should look something like:

 

gnl|Amel_4.0|Group10.1  90      A       10      ........,,      <>[email protected]=BBC      11      ,,,.,.,.,C^~,   CBBCCBBCA-B
gnl|Amel_4.0|Group10.1  91      T       10      ........,,      AABBBCACBB      11      ,,,.,.,.,.,     [email protected]>@=
gnl|Amel_4.0|Group10.1  92      C       10      TTTTTTTTtt      [email protected]?BBB      11      ,,,.,.,.,.,     BBBACA=CB?B
gnl|Amel_4.0|Group10.1  93      G       11      ........,,^~.   [email protected]     11      aaaAaAaAaAa     [email protected]
gnl|Amel_4.0|Group10.1  94      A       11      ........,,.     >;5
gnl|Amel_4.0|Group10.1  96      C       11      TTTTTTTTttT     BAA?BC?CBC;     11      ,,,.,.,.,.,     BBCAABA=BBA
gnl|Amel_4.0|Group10.1  97      A       11      ........,,.     >[email protected]=BCC:     11      ,,,.,.,.,.,     [email protected]
gnl|Amel_4.0|Group10.1  98      G       11      ........,,.     [email protected]<     11      aaaAaAaAaAa     ABBCBABB

Methylation Status report (BED format) for Homogeneous Samples

If the tissue sample being sequenced is expected to be homogenous then you may want to call Methylation state rather than percentage methylated. This mode will report a Cytosine as either Methylated or Unmethylated using the “6 base” model of A,Cm,Cu,Gm,Gu & T. Cytosines are classified as either Methylated or Unmethyated with no shades of grey other than in Diploid mode where it’s possible that a Cytosine that is close to 50 methylated may be called as heterozygous for Cm & Cu.

Here we use novomethyl to produce a bed file showing methylation status of all cytosine’s that had sufficient read coverage to call the status.

novomethyl -h <biseq.mpileup.txt >Cmethyl.bed

The bed file has separate tracks for each context and methylation status. The score column indicates the quality of the call and strand indicates Watson/Crick location of the Cytosine.

gnl|Amel_4.0|Group10.1  91      92      CpGu    54      + 
gnl|Amel_4.0|Group10.1  92      93      CpGu    54      - 
gnl|Amel_4.0|Group10.1  94      95      CHHu    54      + 
gnl|Amel_4.0|Group10.1  95      96      CHGu    54      + 
gnl|Amel_4.0|Group10.1  97      98      CHGu    54      - 
gnl|Amel_4.0|Group10.1  98      99      CpGu    54      + 
gnl|Amel_4.0|Group10.1  99      100     CpGu(Het)    54      - 
gnl|Amel_4.0|Group10.1  101     102     CpGu    54      + 
gnl|Amel_4.0|Group10.1  102     103     CpGu    54      - 
gnl|Amel_4.0|Group10.1  104     105     CHHu    57      - 
gnl|Amel_4.0|Group10.1  106     107     CHHu    57      - 
gnl|Amel_4.0|Group10.1  108     109     CHHu    57      -

A log file is written to stderr with some statistics on Cytosine classifications. The counts of methylated/unmethylated cytosines only includes those where methylation status quality is greater than 5

#       Novomethyl (Beta.4.0) Run Report 
# (C) 2011 NovoCraft Technologies Sdn Bhd. 
# Licensed to Novocraft Technologies Internal Use Only 
# novomethyl -h mpileup.10.3.1M.txt  # Starting at Fri Mar 18 10:13:39 2011 
# 
#       Cytosine Classifications
# 
#       CT strand 
#              660      SNP (C in reference but not in reads.) 
#                0      Low Read Cover
#            14405      Low Quality 
#                0      No Read Cover 
# 
#       Homozygous Cytosines 
#       Context Methylated      Unmethylated    % Methylated 
#       CpG             49           65067             0.1% 
#       CHG              1           24732             0.0% 
#       CHH             12          102977             0.0% 
# 
#       Heterozygous Cytosines #       Context Methylated      Unmethylated    % Methylated 
#       CpG             33             389             7.8% 
#       CHG              1             170             0.6% 
#       CHH              5             646             0.8% 
# 
#       GA strand 
#              536      SNP (C in reference but not in reads.)
#                0      Low Read Cover 
#            13995      Low Quality 
#                0      No Read Cover 
# 
#       Homozygous Cytosines 
#       Context Methylated      Unmethylated    % Methylated 
#       CpG             57           65276             0.1% 
#       CHG              2           24750             0.0% 
#       CHH              5          101362             0.0% 
# 
#       Heterozygous Cytosines 
#       Context Methylated      Unmethylated    % Methylated 
#       CpG             21             359             5.5% 
#       CHG              0             162             0.0% 
#       CHH              1             629             0.2% 
# 
#       99.9%   Cu Conversion Efficiency (CT strand) 
#       99.9%   Cu Conversion Efficiency (GA strand) 
#        Note. Conversion efficiency may be inaccurate if the sample is heterogeneous. 
# Finished at Fri Mar 18 10:13:46 2011

Percent Methylation Report (BED format) for Heterogeneous Samples

Novomethyl is used to produce a bed file showing methylation percentage for each cytosine that had sufficient read coverage to call the status.

 

novomethyl -h -% <biseq.mpileup.txt >Cmethyl.bed

The bed file has separate tracks for each context and methylation status. The score column indicates the % of methylated cytosines and the strand indicates Watson/Crick location of the Cytosine. !

gnl|Amel_4.0|Group10.3  350     350     CHH     0       +
gnl|Amel_4.0|Group10.3  362     362     CHG     0       +
gnl|Amel_4.0|Group10.3  398     398     CHG     0       +
gnl|Amel_4.0|Group10.3  410     410     CHG     40      +
gnl|Amel_4.0|Group10.3  422     422     CHH     0       +
gnl|Amel_4.0|Group10.3  446     446     CHH     0       +
gnl|Amel_4.0|Group10.3  448     448     CHG     0       +

A log file is written to stderr. This includes a histogram of methylation rates for cytosines called with a quality >greater than -q setting and methylation status quality greater then 5.

#       Novomethyl (Beta.4.0) Run Report
# (C) 2011 NovoCraft Technologies Sdn Bhd.
# Licensed to Novocraft Technologies Internal Use Only
# novomethyl -% mpileup.10.3.1M.txt 
# Starting at Fri Mar 18 10:14:02 2011
#
#       Cytosine Classifications
#
#       CT strand
#              210      SNP (C in reference but not in reads.)
#            14799      Low Read Cover
#              264      Low Quality
#                0      No Read Cover
#       % Meth  CpG     CHG     CHH
#       0       65038   24759   103078
#       1       0       0       0
#       2       0       0       0
#       3       0       0       0
#       4       0       0       1
#       5       4       2       0
#       6       18      8       11
.
.
.
#       96      0       0       0
#       97      0       0       0
#       98      0       0       0
#       99      0       0       0
#       100     39      1       10
#
#       GA strand
#              208      SNP (C in reference but not in reads.)
#            14216      Low Read Cover
#              270      Low Quality
#                0      No Read Cover
#       % Meth  CpG     CHG     CHH
#       0       65284   24771   101513
#       1       1       0       0
#       2       0       0       0
#       3       0       0       0
.
.
.
#       97      0       0       0
#       98      0       0       0
#       99      0       0       0
#       100     50      2       5
#
#       99.9%   Cu Conversion Efficiency (CT strand)
#       99.9%   Cu Conversion Efficiency (GA strand)
#        Note. Conversion efficiency may be inaccurate if the sample is heterogeneous.
# Finished at Fri Mar 18 10:14:07 2011

The Cu Conversion efficiencies are tallied from Cytosines that are called as unmethylated.

Consensus & Methylation Report

In consensus report format the novomethyl report is similar to a samtools pileup report with added fields for cytosine context, status, quality and percent methylated.

novomethyl -h -o Consensus< <biseq.mpileup.txt >Cmethyl.txt

An example report is shown below followed by a description of the columns.

gnl|Amel_4.0|Group10.3  1401    T       T       65      .       .       .       .       .       .       7       T$T$TTTTT       =??ABBA 2       TT      CA
gnl|Amel_4.0|Group10.3  1402    T       T       58      .       .       .       .       .       .       5       TTTTT   [email protected]   2       TT      CB
gnl|Amel_4.0|Group10.3  1403    A       A       64      .       .       .       .       .       .       5       AAAAA   @?BBA   2       AA      @A
gnl|Amel_4.0|Group10.3  1404    C       C       56      CHGu    +       66      0       0       5       5       tTttt   BBCAB   2       CC      [email protected]
gnl|Amel_4.0|Group10.3  1405    C       C       56      CpGu    +       66      0       0       5       5       tTttt   [email protected]   2       CC      [email protected]
gnl|Amel_4.0|Group10.3  1406    G       G       63      CpGu    -       55      0       0       2       5       GGGGG   BABCB   2       AA      [email protected]
gnl|Amel_4.0|Group10.3  1407    T       T       58      .       .       .       .       .       .       5       TTTTT   B:ACC   2       TT      BA
gnl|Amel_4.0|Group10.3  1408    A       A       64      .       .       .       .       .       .       5       AAAAA   @[email protected]@   2       AA      [email protected]
gnl|Amel_4.0|Group10.3  1409    T       T       58      .       .       .       .       .       .       5       TTTTT   @ABAB   2       TT      CB
...
gnl|Amel_4.0|Group10.3  1553    T       C/T     66      CHGu    +       27      0       0       4       9       TTTTTTTTT       [email protected]       10      ccTcCcTTCT      [email protected]?
...
gnl|Amel_4.0|Group10.3  88622   A       A       82      .       .       .       .       .       .       8       AAAAAAA^~A      [email protected]        2       AA      A?
gnl|Amel_4.0|Group10.3  88623   A       A       82      .       .       .       .       .       .       8       AAAAAAAA        8;BABBCC        2       AA      C7
gnl|Amel_4.0|Group10.3  88624   C       C       56      CHGu    +       84      0       0       8       8       tTttTTTt        7>[email protected]@C=        2       Cn      B!
gnl|Amel_4.0|Group10.3  88625   C       C       15      CpGm    +       10      75      6       8       8       CCCCTCTC        B2>[email protected]        2       CC      ;>
gnl|Amel_4.0|Group10.3  88626   G       G       84      CpGm    -       32      100     2       2       8       GGGGGGGG        [email protected];        2       GG      AB
gnl|Amel_4.0|Group10.3  88627   A       A       88      .       .       .       .       .       .       9       AAAAAAAA^~A     @B<[email protected]@B       2       AA      B=
gnl|Amel_4.0|Group10.3  88628   C       C       9       CpGu    +       4       67      6       9       9       CCCCTCTtC       B7CBAABCC       2       CC      B;
gnl|Amel_4.0|Group10.3  88629   G       G       90      CpGm    -       32      100     2       2       9       GGGGGGGGG       BAA=ABCBB       2       GG      AA
gnl|Amel_4.0|Group10.3  88630   A       A       88      .       .       .       .       .       .       9       AAAAAAAAA       [email protected]@BAB       2       AA      @?
gnl|Amel_4.0|Group10.3  88631   A       A       88      .       .       .       .       .       .       9       A$AAAAAAAA      A)BB<BBCC       2       AA      BA

A log file is also written to stderr. The log file is a combination of the status and methylation percentage reports. For sake of brevity an example is not given.

 

Consensus Report Columns

1. RNAME Reference sequence name from pileup
2. POS 1-based Coordinate
3. REF Reference Sequence Base. May include IUB ambiguous codes if present in the mpileup.
4. CONSENSUS Consensus base, may show heterozygous base such as A/T
5. QUAL Quality of the consensus call
6. CONTEXT If it is a Cytosine we show the context as CpG, CHG or CHH together with methylation status indicator ‘u’ or ‘m’ if methylation quality is gretaer than 3.
7. STRAND Watson-Crick location of the cytosine is shown as +/-
8. MQUAL The quality of the methylation status given that we have a cytosine at this locus
9. MPCNT The % of cytosines that were not converted to T’s or a period if the number of samples is too low to accurately call a percentage.
10. NU The number of unconverted cytosines.
11. NC The total number of cytosines. If a site is called as heterozygous C/T then the lesser of half the total C&T bases or the total T bases are assumed to result from the heterozygous T and are not included in the NC count.
12. NCT The number of reads aligning on the Watson strand. The following 6 fields are directly from samtools mpileup file and we refer you to samtools documentation for further information.
13. CTBASES The bases and read mapping qualities from Watson strand.
14. CTQUALITIES The base qualities for Watson strand.
15. NGA Next 3 fields are as above for Crick strand.
16. GABASES
17. GAQUALITIES

Novomethyl Command Line Option

Usage:

 

novomethyl options < mpileup.txt Options (with default value shown):

-1 Sets Haploid mode. The sample is expected to be from a Haploid genome and and heterozygosity is not considered as a possible base state.
-2 Sets Diploid mode (default) and only reports Homozygous Cytosines. The sample is expected to be from a Diploid genome and and heterozygosity is considered as a possible base state. In BED format, only Homozygous cytosines have their methylation status or percentage reported.
-h Turns on Diploid mode with reporting of both homozygous and heterozygous Cytosines. The sample is expected to be from a Diploid genome and and heterozygosity is considered as a possible base state. Homozygous & Heterozygous cytosines will have their methylation status or percentage reported.
-% In BED format report % methylated rather than methylation status and quality. In this mode we are calling % methylation for the cytosine. This can be combined with the above Haploid & Diploid modes. If in Diploid mode (-h) and a locii is identified as Heterozygous C/[AGT] then % methylation is calculated on the C  and can still show 100% methylated.
-o Consensus Uses consensus sequence report format. This format contains much more information however you may have to filter and process it to make a track suitable for your genome browser.
-s 0.001 Sets prior probability of a SNP. Used to set prior probabilities that a base is different to the reference genome. The reference base gets a prior of 1.0 all other bases get prior of Psnp. Priors also take into account ambiguous nucleotides in the reference genome.
-e 0.98 Sets bisulphite conversion efficiency. Used in calculation of P(dk|H) for unconverted Cytosines. Refer to SOAP paper for more details. We suggest running a test and then adjusting this to match the reported rates especially if your sample is heterogeneous.
-m 0.5 Sets prior probability of Cytosine being methylated. Splits prior probability of Cytosine between methylated and unmethylated states so that prior of a methylated Cytosine is P(C).Pm
-q 10 Sets quality limit for reporting Cytosine methylation status.(Note. A quality of 3 is probability of 0.5 that Cytosine is methylated so it’s really an unknown status.) In both state and % mode a cytosine locii will not be reported if it’s quality is less than this threshold. In state mode this is quality of cm or Cu state, in percent mode it is the quality of the C call. Quality is calculated using a 6 base model (A,Cu,Cm,Gu,Gm,T) with the basic statistical models from SOAP SNP detection for massively parallel whole-genome resequencing” (2009) Genome Res. published online May 6, 2009 , doi:10.1101/gr.088013.108
-Q 10 Lower limit on base quality for counting bases for % methylation and conversion rates. Base calls from  mpileup file with quality less than this threshold will not be counted for determining % methylation or the conversion efficiency.
-c 1 Sets minimum number of reads covering a Cytosine (strand specific). In percent mode a Cytosine will not be called unless it has at least this minimum number of reads on CT or GA strand.
-d dbname Indexed reference genome(optional). This is the novoindex of the reference genome used to produce the alignments. If provided it is used for determing the cytosine context (CpG, CHG, CHH).It helps resolve context near the edges of read cover where the mpileup file does not include all the references base need to establish the context.
-D Debug mode (prints pileup details) The Consensus format report is possibly more useful.

Heterozygous C/G Locii

If a locus is identified as heterozygous for C&G Novomethyl will report the site only once. The consensus is identified as C/G with appropriate quality and then methylation state and percentage is reported for either the C or the G but not both. The strand column will indicate if it is the Watson or Crick cytosine that is reported. The strand that has the higher quality for the cytosine (C or G) call is reported.

 

default
Loading posts...
link_magnifier
#6E787E
on
fadeInDown
loading
#6E787E
on