Error estimation - Revision history

Thorfinn: /* Error estimation using an outgroup and an error free individual */

2017-01-12T11:57:39Z

Error estimation using an outgroup and an error free individual

← Older revision		Revision as of 13:57, 12 January 2017
Line 90:		Line 90:

	The estimated rates can roughly be interpreted as relative error rates. That is excess of derived alleles in your sample compare to the derived alleles in the an error free indviduals The idea is the your sample and the error free individuals should have the same expected number of derived alleles and the extra observed derived alleles in you sample are due to the excess error. For each individual we sample a single base from the reads at each position. We use only positions were there are coverage for both the outgroup, the sample and the error free individual. Some more details about the method is described [[ErrorPerfect\|here]].		The estimated rates can roughly be interpreted as relative error rates. That is excess of derived alleles in your sample compare to the derived alleles in the an error free indviduals The idea is the your sample and the error free individuals should have the same expected number of derived alleles and the extra observed derived alleles in you sample are due to the excess error. For each individual we sample a single base from the reads at each position. We use only positions were there are coverage for both the outgroup, the sample and the error free individual. Some more details about the method is described [[ErrorPerfect\|here]].

			If you work with human data mapped to hg19 we have a fasta file containing the ancestral states (actually the chimp mapped to hg19).
			<pre>
			wget http://popgen.dk/software/download/angsd/hg19ancNoChr.fa.gz
			</pre>
	==Brief Overview==		==Brief Overview==
	<pre>		<pre>

Thorfinn: /* Output */

2015-11-26T12:22:06Z

Output

← Older revision		Revision as of 14:22, 26 November 2015
Line 153:		Line 153:
	==Output==		==Output==
	# The output from the angsd cprogram is in a horrible format. Those are 2 files suffixed with '''.ancError''' and '''.ancErrorChr'''.		# The output from the angsd cprogram is in a horrible format. Those are 2 files suffixed with '''.ancError''' and '''.ancErrorChr'''.
			For each individual we allocate a integer vector with 125 entries. Each entry represents the counts of combinations of the outgroup,perfect,sample alleles. The offset is given by: outgroup25+perfect5+bam;With the bases encoded as a=0,c=1,g=2,t=3,n=4

	# The Rscript will output the type specific errors and the overall errors. It will also produce two figures. One will the type specific errors and one with the overall errors. If a lot of individuals are included in the analysis the figure will need to be modified. The output is 3 files: 2) figures called '''errorEst/errorEstOverall'''. And a text file containing the values used for generating the figures. Example below		# The Rscript will output the type specific errors and the overall errors. It will also produce two figures. One will the type specific errors and one with the overall errors. If a lot of individuals are included in the analysis the figure will need to be modified. The output is 3 files: 2) figures called '''errorEst/errorEstOverall'''. And a text file containing the values used for generating the figures. Example below

Line 168:		Line 170:
	"ind5 0.019 %"		"ind5 0.019 %"
	</pre>		</pre>


	=Mismatch rates=		=Mismatch rates=

	Type, strand, read position specific mismatch rates can be obtain using the [[Genotype_likelihoods#SOAPsnp \| soapSNP genotype likelihood method]]. Use a reference genome with known variations removed to reduce the bias of this method		Type, strand, read position specific mismatch rates can be obtain using the [[Genotype_likelihoods#SOAPsnp \| soapSNP genotype likelihood method]]. Use a reference genome with known variations removed to reduce the bias of this method

Albrecht at 15:25, 26 February 2014

2014-02-26T15:25:34Z

← Older revision		Revision as of 17:25, 26 February 2014
Line 168:		Line 168:
	"ind5 0.019 %"		"ind5 0.019 %"
	</pre>		</pre>


			=Mismatch rates=

			Type, strand, read position specific mismatch rates can be obtain using the [[Genotype_likelihoods#SOAPsnp \| soapSNP genotype likelihood method]]. Use a reference genome with known variations removed to reduce the bias of this method

Albrecht: /* Theory */

2014-02-26T15:23:17Z

Theory

← Older revision		Revision as of 17:23, 26 February 2014
Line 168:		Line 168:
	"ind5 0.019 %"		"ind5 0.019 %"
	</pre>		</pre>

	~~=Theory=~~
	~~These two subsection will contain the theory behind the methods.~~

Albrecht: /* Example */

2014-02-26T15:21:29Z

Example

← Older revision		Revision as of 17:21, 26 February 2014
Line 146:		Line 146:

	<pre>		<pre>
	Rscript ~~../angsd0.567/~~R/estError.R file=output/smallBamMicsAncq30Q20.ancError main="Errorrate using outgroup (CHIMP) and a perfect MAN" doPng=TRUE		Rscript R/estError.R file=output/smallBamMicsAncq30Q20.ancError main="Errorrate using outgroup (CHIMP) and a perfect MAN" doPng=TRUE
	</pre>		</pre>

	;NB the script will generate an empty Rplots.pdf file. This is expected. Sorry about this. The reason is that we haven't figured out how to change the R command called palette without generating a Rplots.pdf file.		;NB the script will generate an empty Rplots.pdf file. This is expected. Sorry about this. The reason is that we haven't figured out how to change the R command called palette without generating a Rplots.pdf file (stupid R 3.+ color palette).

	==Output==		==Output==

Albrecht: /* Options */

2014-02-26T15:19:10Z

Options

← Older revision		Revision as of 17:19, 26 February 2014
Line 117:		Line 117:

	==Options==		==Options==
	; -doAncError 1		; -doAncError [int]
	Use -doAncError to use this analysis		Use -doAncError to use this analysis.
	~~; -anc [filename]~~
	~~fasta file with the ancestral alleles~~
	~~; -ref [filename]~~
	~~fasta file of an error free individual~~.
	~~; -sample [int]~~
	0: No error estimation		0: No error estimation

Line 130:		Line 125:
	2: Sample only one base per site		2: Sample only one base per site

	3: Use first ~~base~~ per site		3: Use first read per site (the reads are order with position start)

	To further refine what data should be used please see [[alleles counts]] and [[Input \| read filters]].		To further refine what data should be used please see [[alleles counts]] and [[Input \| read filters]].


			; -anc [filename]
			fasta file with the ancestral alleles
			; -ref [filename]
			fasta file of an error free individual.
			]

	==Example==		==Example==

Albrecht: /* Options */

2014-02-26T15:18:02Z

Options

← Older revision		Revision as of 17:18, 26 February 2014
Line 124:		Line 124:
	fasta file of an error free individual.		fasta file of an error free individual.
	; -sample [int]		; -sample [int]
	0: ~~Use all reads~~		0: No error estimation

	1: ~~Sample only one base per site (default)~~		1: Use all reads

	2: Use first base per site		2: Sample only one base per site

			3: Use first base per site

	To further refine what data should be used please see [[alleles counts]] and [[Input \| read filters]].		To further refine what data should be used please see [[alleles counts]] and [[Input \| read filters]].

Albrecht: /* Error estimation using an outgroup and an error free individual */

2014-02-26T15:16:44Z

Error estimation using an outgroup and an error free individual

← Older revision		Revision as of 17:16, 26 February 2014
Line 93:		Line 93:
	<pre>		<pre>
	./angsd -doAncError		./angsd -doAncError
	-> angsd version: 0.~~567~~ build(~~Dec 7 2013 14~~:56:25)		Command:
			./angsd -doAncError
			-> angsd version: 0.580 build(Feb 26 2014 10:11:49)
	-> Analysis helpbox/synopsis information:		-> Analysis helpbox/synopsis information:
	--------------		--------------
	analysisAncError.cpp:		analysisAncError.cpp:
	-doAncError 0		-doAncError 0
	~~1: EM v1~~		(Sampling strategies)
	~~2: EM v2 (Under development, don't use)~~		0: no error estimation
	~~-sample 1~~ (Sampling strategies~~, for EM v1~~)		1: (Use all bases)
	0: (Use all bases)		2: (Sample single base)
	1: (Sample single base)		3: (Sample first base)
	2: (Sample first base)
	Required:		Required:
	-ref (null) (fastafile containg 'perfect' sample)		-ref (null) (fastafile containg 'perfect' sample)
	-anc (null) (fastafile containg outgroup)		-anc (null) (fastafile containg outgroup)

			NB: the -ref should be a fasta for a sample where you assume no errors.
			We measure the difference between the outgroup and your -ref sample.
			The statistic is then the excess of substitutions between your BAM file and outgroup, compared to the perfect sample. After the ANGSD run use: Rscript R/estError.R file=angsdput.ancerror
	</pre>		</pre>

Albrecht: /* Error estimation using an outgroup and an error free individual */

2014-02-26T15:13:56Z

Error estimation using an outgroup and an error free individual

← Older revision		Revision as of 17:13, 26 February 2014
Line 89:		Line 89:
	[[File:errorEstOverall.png\|thumb]]		[[File:errorEstOverall.png\|thumb]]

	The estimated rates can roughly be interpreted as relative error rates. That is excess of derived alleles in your sample compare to the derived alleles in the an error free indviduals The idea is the your sample and the error free individuals should have the same expected number of derived alleles and the extra observed derived alleles in you sample are due to the excess error. For each individual we sample a single base from the reads at each position. We use only positions were there are coverage for both the outgroup, the sample and the error free individual. ~~See theory section below for~~ more details.		The estimated rates can roughly be interpreted as relative error rates. That is excess of derived alleles in your sample compare to the derived alleles in the an error free indviduals The idea is the your sample and the error free individuals should have the same expected number of derived alleles and the extra observed derived alleles in you sample are due to the excess error. For each individual we sample a single base from the reads at each position. We use only positions were there are coverage for both the outgroup, the sample and the error free individual. Some more details about the method is described [[ErrorPerfect\|here]].
	==Brief Overview==		==Brief Overview==
	<pre>		<pre>

Albrecht: /* Method 2 */

2014-02-26T15:13:11Z

Method 2

← Older revision		Revision as of 17:13, 26 February 2014
Line 161:		Line 161:

	=Theory=		=Theory=
	These two subsection will contain the theory behind the methods.		These two subsection will contain the theory behind the methods.
	~~==Method 2==~~

	~~The method was use in [[Orlando2013]] and is described in details in supplementary 4.4~~

	The estimated rates can roughly be intrepreted as relative error rates. That is excess of errors in your sample compare to the error in a high quality indvidual. The idea is the your sample and the high quality individual should have the same expected number of derived alleles and the extra derived alleles in you sample are due to the excess error. We use only positions were there are coverage for both the chimp, the sample and the high quality genome. The overall error rate is obtained from


	~~<math>~~
	~~O_D = E_D(1-\epsilon) + E_A\epsilon~~
	~~</math>~~

	~~were~~

	*<math>\epsilon</math> is the error rate
	*<math>O_D</math> is the observed number of derived alleles in the sample
	*<math>E_D</math> is the expected number of derived alleles which is obtained from the observed derived alleles from the high quality genome
	*<math>E_A</math> is the expected number of ancestral alleles which is obtained from the high quality genome




	~~The type specific error rates are obtained from maximizing the likehood~~

	~~<math>~~
	~~P(H=h\|D=d,e) = \prod_1^n P(H_i=h_i\|D_i=d_i,e)~~
	~~</math>~~
	where n is the number of sampled bases, <math>h_i </math> is the observed base of the sample, <math> d_i</math> is the obsered base of the outgroup and e is an error rate matrix with entrance <math>e_{a->b} </math> being the rate of
	~~error from base a to different base b and <math>e_{a->a} </math> being equal to <math>1-\sum_{b\ne a} e_{a->b}</math>~~

	~~We calculate the above by summing over the true unobserved base <math>T_{i} </math>~~
	~~<math>\begin{align}~~
	~~P(H_i=h_i\|D_i=d_i,e) &= \sum_{t_i\in\{A,C,G,T,\}} P(H_i=h_i,T_i=t_i\|D_i=d_i,e) \\~~
	~~&= \sum_{t_i\in\{A,C,G,T,\}} P(H_i=h_i\|T_i=t_i,e)P(T_i=t_i\|D_i=d_i) \\~~
	~~&= \sum_{t_i\in\{A,C,G,T,\}} e_{t_i->h_i}P(T_i=t_i\|D_i=d_i)~~
	~~\end{align}~~
	~~</math>~~


	<math>P(T_i=t_i\|D_i=d_i) </math> is set equal to the observed fraction of sites in the high quality genome where the outgroup has base <math>d_i</math> and the high qualtiy genome has base <math>t_i</math>