Safv3 - Revision history

Thorfinn at 15:47, 5 February 2016

2016-02-05T15:47:25Z

← Older revision		Revision as of 17:47, 5 February 2016
Line 2:		Line 2:

	This page will contain the impact of this new format in downstream analysis.		This page will contain the impact of this new format in downstream analysis.

			NB the ml estimate from realSFS is now expectations of the sfs, and not in logscale. So the exp commands in the code belove is deprecated.

Thorfinn: /* One population analysis */

2016-02-05T15:46:05Z

One population analysis

← Older revision		Revision as of 17:46, 5 February 2016
Line 71:		Line 71:
	<pre>		<pre>
	norm<-function(X) X/sum(X)		norm<-function(X) X/sum(X)
	a<-~~exp(~~read.table("oldceu2.oldyri2.100mb",nrow=37))		a<-read.table("oldceu2.oldyri2.100mb",nrow=37)
	b<-norm(scan("newceu.newyri.chr1.100mb"))		b<-norm(scan("newceu.newyri.chr1.100mb"))
	b<-matrix(b,ncol=53,byrow=T)		b<-matrix(b,ncol=53,byrow=T)
Line 139:		Line 139:
	}		}
	if(FALSE){		if(FALSE){
	pop1 <- ~~exp(~~scan("pop1.saf.idx.ml"))		pop1 <- scan("pop1.saf.idx.ml")
	pop2 <- ~~exp(~~scan("pop2.saf.idx.ml"))		pop2 <- scan("pop2.saf.idx.ml")
	pop1.pop2 <- matrix(exp(scan("pop1.pop2.saf.idx.ml")),nPop1+1,byrow=T)		pop1.pop2 <- matrix(exp(scan("pop1.pop2.saf.idx.ml")),nPop1+1,byrow=T)
	par(mfrow=c(1,2))		par(mfrow=c(1,2))

Thorfinn: /* One population analysis */

2016-02-05T15:45:16Z

One population analysis

← Older revision		Revision as of 17:45, 5 February 2016
Line 21:		Line 21:

	##comparison		##comparison
	a<-~~exp(~~scan("newceu.saf.idx.chr1.ml"))		a<-scan("newceu.saf.idx.chr1.ml")
	b<-~~exp(~~as.numeric(read.table("oldceu.saf.ml")[1,]))		b<-as.numeric(read.table("oldceu.saf.ml")[1,])
	a-b		a-b
	[1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00		[1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

Thorfinn: /* Two population analysis */

2015-05-11T08:09:19Z

Two population analysis

← Older revision		Revision as of 10:09, 11 May 2015
Line 76:		Line 76:
	> range(a-b)		> range(a-b)
	[1] -1.724603e-07 2.850303e-09		[1] -1.724603e-07 2.850303e-09

			> png("safv3to.png")
			> barplot(rbind(colSums(a),colSums(b))[,-c(1,53)],be=T,main="marginal 2dsfs YRI",legend=c("safv1","safv3"),col=1:2)
			> dev.off()
			X11cairo
			2
			> png("safv3en.png")
			> barplot(rbind(rowSums(a),rowSums(b))[,-c(1,37)],be=T,main="marginal 2dsfs CEU",legend=c("safv1","safv3"),col=1:2)
			> dev.off()

	</pre>		</pre>
	[[File:safv3en.png\|thumb]]		[[File:safv3en.png\|thumb]]
	[[File:safv3to.png\|thumb]]		[[File:safv3to.png\|thumb]]
			This looks very nice, right?

	==2 population analysis with simulated data==		==2 population analysis with simulated data==

Thorfinn: /* Two population analysis */

2015-05-11T07:55:53Z

Two population analysis

← Older revision		Revision as of 09:55, 11 May 2015
Line 77:		Line 77:
	[1] -1.724603e-07 2.850303e-09		[1] -1.724603e-07 2.850303e-09
	</pre>		</pre>
			[[File:safv3en.png\|thumb]]
			[[File:safv3to.png\|thumb]]

	==2 population analysis with simulated data==		==2 population analysis with simulated data==

Thorfinn: /* Two population analysis */

2015-05-11T07:50:30Z

Two population analysis

← Older revision		Revision as of 09:50, 11 May 2015
Line 70:		Line 70:
	and comparing in R		and comparing in R
	<pre>		<pre>
			norm<-function(X) X/sum(X)
			a<-exp(read.table("oldceu2.oldyri2.100mb",nrow=37))
			b<-norm(scan("newceu.newyri.chr1.100mb"))
			b<-matrix(b,ncol=53,byrow=T)
			> range(a-b)
			[1] -1.724603e-07 2.850303e-09
	</pre>		</pre>

Thorfinn: /* Two population analysis */

2015-05-11T07:46:39Z

Two population analysis

← Older revision		Revision as of 09:46, 11 May 2015
Line 67:		Line 67:

	../angsd/misc/realSFS ../nsfs/newceu.saf.idx ../nsfs/newyri.saf.idx -P 32 -r 1 -nSites 100000000 >ceu.yri.chr1.100mb		../angsd/misc/realSFS ../nsfs/newceu.saf.idx ../nsfs/newyri.saf.idx -P 32 -r 1 -nSites 100000000 >ceu.yri.chr1.100mb
			</pre>
			and comparing in R
			<pre>

	</pre>		</pre>

Thorfinn: /* Two population analysis */

2015-05-11T07:44:20Z

Two population analysis

← Older revision		Revision as of 09:44, 11 May 2015
Line 45:		Line 45:
	==> oldceu2.arg <==		==> oldceu2.arg <==
	../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu2 -r 1 -sites intersect.txt		../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu2 -r 1 -sites intersect.txt

	==> oldceu.arg <==		==> oldceu.arg <==
	../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu -~~rf rf~~		../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu -r 1

	==> oldyri2.arg <==		==> oldyri2.arg <==
	../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri2 -r 1 -sites intersect.txt		../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri2 -r 1 -sites intersect.txt
			==> oldyri.arg <==
			../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri -r 1

	~~==> oldyri.arg <==~~
	~~../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri -rf rf~~
	~~==> oldceu2.arg <==~~
	~~../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu2 -rf rf -sites intersect.txt~~
	~~==> oldceu.arg <==~~
	~~../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu -rf rf~~
	~~==> oldyri2.arg <==~~
	~~../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri2 -rf rf -sites intersect.txt~~
	~~==> oldyri.arg <==~~
	~~../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri -rf rf~~

	##with intersect found like		##with intersect found like
Line 67:		Line 57:


	##The old saf files are very big so we had to limit the analysis to ~~180mio~~ sites		##The old saf files are very big so we had to limit the analysis to 100mio sites
	../master/misc/realSFS 2dsfs oldceu2.saf oldyri2.saf 36 52 -nSites ~~180000000~~ -P 20 >oldceu2.oldyri2.ml		../master/misc/realSFS 2dsfs oldceu2.saf oldyri2.saf 36 52 -nSites 100000000 -P 20 >oldceu2.oldyri2.ml.100mb

	##the new format is much simpler here we simply did		##the new format is much simpler here we simply did
	==> newceu.arg <==		==> newceu.arg <==
	../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out newceu -~~rf rf~~		../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out newceu -r 1
	==> newyri.arg <==		==> newyri.arg <==
	../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out newyri -~~rf rf~~		../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out newyri -r 1

	../angsd/misc/realSFS ../nsfs/newceu.saf.idx ../nsfs/newyri.saf.idx -P 32 -r 1 >ceu.yri.chr1.~~180mb~~		../angsd/misc/realSFS ../nsfs/newceu.saf.idx ../nsfs/newyri.saf.idx -P 32 -r 1 -nSites 100000000 >ceu.yri.chr1.100mb
	</pre>		</pre>

Thorfinn: /* Two population analysis */

2015-05-11T07:41:39Z

Two population analysis

← Older revision		Revision as of 09:41, 11 May 2015
Line 39:		Line 39:
	==Two population analysis==		==Two population analysis==
	<pre>		<pre>
			##old(master) -> angsd version: 0.801-28-gbab908a (htslib: 1.2.1-62-g35746af) build(May 9 2015 14:50:33)
			##new(newsaf)-> angsd version: 0.801-61-g48f06d8-dirty (htslib: 1.2.1-62-g35746af) build(May 9 2015 08:49:40)
	##old version required a run for each population, to find the intersect and then limit the analysis to the intersect.		##old version required a run for each population, to find the intersect and then limit the analysis to the intersect.
	##Here are all 4 commands		##Here are all 4 commands
			==> oldceu2.arg <==
			../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu2 -r 1 -sites intersect.txt

			==> oldceu.arg <==
			../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu -rf rf

			==> oldyri2.arg <==
			../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri2 -r 1 -sites intersect.txt

			==> oldyri.arg <==
			../angsd/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/yri.ricco.list -gl 1 -P 5 -out oldyri -rf rf
	==> oldceu2.arg <==		==> oldceu2.arg <==
	../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu2 -rf rf -sites intersect.txt		../master/angsd -anc hg19ancNoChr.fa -dosaf 1 -b /space/genomes/1000g/lowC2014/filelists/ceu.ricco.list -gl 1 -P 5 -out oldceu2 -rf rf -sites intersect.txt

Thorfinn: /* 2 population analysis with simulated data */

2015-05-07T13:34:43Z

2 population analysis with simulated data

← Older revision		Revision as of 15:34, 7 May 2015
Line 71:		Line 71:

	<pre>		<pre>
	##~~taken from~~ angsd~~/R/2dsfs~~.R		##angsd new: ->-> angsd version: 0.801-51-g156039a (htslib: 1.2.1-69-gb79f40a) build(May 7 2015 15:31:53)
			##angsd old: -> angsd version: 0.801-27-ga699b44 (htslib: 1.2.1-69-gb79f40a) build(May 7 2015 15:30:06)

	norm <- function(x) x/sum(x)		norm <- function(x) x/sum(x)
	if(FALSE){		if(FALSE){
Line 143:		Line 145:
	pop1.pop2 <- matrix(exp(scan("pop1.pop2.saf.idx.ml")),nPop1+1,byrow=T)		pop1.pop2 <- matrix(exp(scan("pop1.pop2.saf.idx.ml")),nPop1+1,byrow=T)
	par(mfrow=c(1,2))		par(mfrow=c(1,2))
	barplot(rbind(norm(table(p1.d)[-1]),pop1),be=T,main="~~only~~ varsites pop1")		barplot(rbind(norm(table(p1.d)[-1]),pop1),be=T,main="varsites pop1")
	barplot(rbind(norm(table(p2.d)[-1]),pop2),be=T,main="~~only~~ varsites pop2")		barplot(rbind(norm(table(p2.d)[-1]),pop2),be=T,main="varsites pop2")
			pop1.pop2[1,1] <- 0
			pop1.pop2[nrow(pop1.pop2),ncol(pop1.pop2)] <- 0
			pop1.pop2 <- norm(pop1.pop2)
	range(norm(sfs.2d)-t(pop1.pop2))		range(norm(sfs.2d)-t(pop1.pop2))
	## [1] -0.0005047007 0.0007018686		#[1] -0.0005047007 0.0007018686

	barplot(rbind(rowSums(norm(t(sfs.2d))),rowSums(pop1.pop2)),be=T)		barplot(rbind(rowSums(norm(t(sfs.2d))),rowSums(pop1.pop2)),be=T)
	barplot(rbind(colSums(norm(t(sfs.2d))),colSums(pop1.pop2)),be=T)		barplot(rbind(colSums(norm(t(sfs.2d))),colSums(pop1.pop2)),be=T)

	}		}
			if(FALSE){
			##just redo the angsd and optimization
			##git checkout master ;make clean;make
			system("../angsd -glf pop1.glf.gz -nind 12 -doSaf 1 -out pop1 -fai fai.fai -issim 1 -P 10")
			system("../angsd -glf pop2.glf.gz -nind 8 -doSaf 1 -out pop2 -fai fai.fai -issim 1 -P 10")
			system("../misc/realSFS pop1.saf 24 -P 60 >pop1.saf.idx.ml")
			system("../misc/realSFS pop2.saf 16 -P 60 >pop2.saf.idx.ml")
			system("../misc/realSFS 2dsfs pop1.saf pop2.saf 24 16 -P 60 >pop1.pop2.saf.idx.ml")
			}
			if(FALSE){
			pop1 <- norm(exp(scan("pop1.saf.idx.ml"))[-1])
			pop2 <- norm(exp(scan("pop2.saf.idx.ml"))[-1])
			pop1.pop2 <- matrix(exp(scan("pop1.pop2.saf.idx.ml")),nPop1+1,byrow=T)
			par(mfrow=c(1,2))
			barplot(rbind(norm(table(p1.d)[-1]),pop1),be=T,main="varsites pop1")
			barplot(rbind(norm(table(p2.d)[-1]),pop2),be=T,main="varsites pop2")
			pop1.pop2[1,1] <- 0
			pop1.pop2[nrow(pop1.pop2),ncol(pop1.pop2)] <- 0
			pop1.pop2 <- norm(pop1.pop2)
			range(norm(sfs.2d)-t(pop1.pop2))
			## [1] -0.0005046856 0.0007019217

			}

	}		}



	</pre>		</pre>