ANGSD: Analysis of next generation Sequencing Data
Latest tar.gz version is (0.938/0.939 on github), see Change_log for changes, and download it here.
Contamination: Difference between revisions
No edit summary |
No edit summary |
||
| Line 11: | Line 11: | ||
<pre> | <pre> | ||
./angsd -i my.bam -r X: -doCounts 1 -iCounts 1 -minMapQ 30 -minQ 20 | ./angsd -i my.bam -r X: -doCounts 1 -iCounts 1 -minMapQ 30 -minQ 20 | ||
Rscript | Rscript contamination.R mapFile="map100.chrX.bz2" hapFile="hapMapCeuXlift.map.bz2" countFile="angsdput.icnts.gz" mc.cores=24 | ||
</pre> | </pre> | ||
The '''contamination.R''' program is found in the '''R/''' subfolder, and the resource files are found in the '''RES''' folder. The jackknive procedure can be quite slow, so we allocate 24 cores for this analysis '''mc.cores=24'''. | |||
=Output= | |||
The output from the above command is shown below | |||
<pre> | <pre> | ||
Revision as of 11:58, 27 June 2014
Angsd can estimate contamination, but only for chromosomes that exists in one genecopy (eg chrX for males). This method requires a list of HapMap sites along with their frequency and we also recommend to discard regions with low mappability.
We have included a mappability and HapMap files for chrX these are found in the RES subfolder of the angsd source package. So if you are working with humans, and your sample is a male then you can estimate the contamination with the follow two commands.
- First we generate a binary count file for chrX for a single BAM file (ANGSD cprogram)
- Then we do a fisher test for finding a p-value, and jackknife to get an estimate of contamination (Rprogram)
An example are found below:
./angsd -i my.bam -r X: -doCounts 1 -iCounts 1 -minMapQ 30 -minQ 20 Rscript contamination.R mapFile="map100.chrX.bz2" hapFile="hapMapCeuXlift.map.bz2" countFile="angsdput.icnts.gz" mc.cores=24
The contamination.R program is found in the R/ subfolder, and the resource files are found in the RES folder. The jackknive procedure can be quite slow, so we allocate 24 cores for this analysis mc.cores=24.
Output
The output from the above command is shown below
Rscript ../R/contamination.R mapFile="../RES/map100.chrX.bz2" hapFile="../RES/hapMapCeuXlift.map.bz2" countFile="/space/anders/ida/idaSjov/kostenkitest/contamination/out/V1countKostinki.USER.bam.X.gz" mc.cores=24
Loading required package: multicore
-----------------------
Doing Fisher exact test for Method1:
[,1] [,2]
[1,] 246 157
[2,] 17700 143407
Fisher's Exact Test for Count Data
data: mat
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
10.34000 15.61672
sample estimates:
odds ratio
12.6959
-----------------------
Doing Fisher exact test for Method2:
[,1] [,2]
[1,] 91 55
[2,] 7355 59513
Fisher's Exact Test for Count Data
data: mat2
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
9.466476 19.085589
sample estimates:
odds ratio
13.38675
----------------------
Running jackknife for Method1 (could be slow)
Running jackknife for Method2 (could be slow)
$est
Method1 Method2
Contamination 0.03837625 0.03380983
llh 1034.078 483.5145
SE 0.002630455 0.003900376
$err
[1] 0.01370779
$c
[1] 0.001093589
$est
Method1 Method2
Contamination 0.03837625 0.03380983
llh 1034.078 483.5145
SE 0.002630455 0.003900376
$err
[1] 0.01370779
$c
[1] 0.001093589