| Previous changeset 0:6b7107812931 (2015-07-02) |
|
Commit message:
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac |
|
added:
MLMM.pl MLMM.sh MLMM.xml source_library/emma.mlmm.r source_library/emma.r source_library/mlmm.r source_library/mlmm1.r source_library/plot_MLMM.Ago.r test-data/PCs.txt test-data/genot.txt test-data/kinship.txt test-data/map.txt test-data/output.txt test-data/phenot.txt test-data/rss.txt test-data/step_table.txt |
|
removed:
mlmm/MLMM.pl mlmm/MLMM.sh mlmm/MLMM.xml mlmm/source_library/emma.mlmm.r mlmm/source_library/emma.r mlmm/source_library/mlmm.r mlmm/source_library/mlmm1.r mlmm/source_library/plot_MLMM.Ago.r |
| b |
| diff -r 6b7107812931 -r 380b364980f9 MLMM.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MLMM.pl Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| @@ -0,0 +1,154 @@ +#!/usr/bin/perl + +use strict; +use Getopt::Long; +use Bio::SeqIO; + + +my $usage = qq~Usage:$0 <args> [<opts>] +where <args> are: + -g, --geno <Genotype input> + -i, --info <SNP information. Genome position.> + -p, --pheno <Phenotype input> + -o, --out <output name> + -d, --directory <directory for MLMM R libraries> + -s, --step_number <number of steps. Maximum: 20. Default: 10> + -m, --method <Method: mbonf or extBIC. Default: mbonf> +~; +$usage .= "\n"; + +my ($geno,$map,$pheno,$out,$dir,$steps,$method); + + +GetOptions( + "geno=s" => \$geno, + "info=s" => \$map, + "pheno=s" => \$pheno, + "out=s" => \$out, + "dir=s" => \$dir, + "steps=s" => \$steps, + "method=s" => \$method +); + + +die $usage + if ( !$geno || !$map || !$pheno || !$out || !$dir || !$steps || !$method); + +my $max_steps = 10; +my $plot_opt = "mbonf"; +if ($method && $method ne 'mbonf' && $method ne 'extBIC') +{ + print "Aborted: Method must be mbonf or extBIC.\n"; + exit; +} +else +{ + $plot_opt = $method; +} +if ($steps && $steps !~/\d+/ && $steps > 20 && $steps < 2) +{ + print "Aborted: Number of steps must be greater than 2 and lower than 20.\n"; + exit; +} +else +{ + $max_steps = $steps; +} + + +my $chunk = 2; + +my $RSCRIPT_EXE = "Rscript"; +my $R_DIR = $dir; + + +my $head_trait = `head -1 $pheno`; +my @headers_traits = split(/\t/,$head_trait); +my $trait_name = $headers_traits[1]; + + +open( my $RCMD, ">rscript" ) or throw Error::Simple($!); + + + + +print $RCMD "Y_file <- \"" . $pheno . "\"\n"; +print $RCMD "X_file <- \"" . $geno . "\"\n"; +if($map) +{ + print $RCMD "map_file <- \"$map\"\n"; + print $RCMD "map <- read.table(map_file, sep = \"\\t\", header = T)\n"; +} +print $RCMD "mlmm_data = list()\n"; +print $RCMD "mlmm_data\$chunk <- $chunk\n"; +print $RCMD "mlmm_data\$maxsteps <- $max_steps\n"; +print $RCMD "genot <- read.table(X_file, sep = \"\\t\", header = T)\n"; +print $RCMD "genot_mat <- as.matrix(genot[, 2:ncol(genot)])\n"; +print $RCMD "rownames(genot_mat) <- genot\$Ind_id\n"; + +print $RCMD "phenot <- read.table(Y_file, sep = \"\\t\", header = T)\n"; + + + +# missing data imputation +print $RCMD "genot_imp <- genot_mat\n"; +print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n"; +print $RCMD "for (i in 1:ncol(genot_imp)){genot_imp[is.na(genot_imp[,i]), i] <- average[i]}\n"; + +# kinship matrix computation +print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n"; +print $RCMD "stdev <- apply(genot_imp, 2, sd)\n"; +print $RCMD "genot_stand <- sweep(sweep(genot_imp, 2, average, \"-\"), 2, stdev, \"/\")\n"; +print $RCMD "K_mat <- (genot_stand %*% t(genot_stand)) / ncol(genot_stand)\n"; +print $RCMD "write.table(K_mat, '$out.kinship', sep='\\t', dec='.', quote=F, col.names=T, row.names=T)\n"; + +print $RCMD "source(\"" . $R_DIR. "/mlmm.r\")\n"; +print $RCMD "source(\"" . $R_DIR. "/emma.r\")\n"; + +# mlmm +print $RCMD "mygwas <- mlmm(Y = phenot\$$trait_name, X = genot_imp, K = K_mat, nbchunks=mlmm_data\$chunk, maxsteps=mlmm_data\$maxsteps)\n"; + +# plots +print $RCMD "pdf('$out.pdf')\n"; +print $RCMD "plot_step_table(mygwas, \"h2\")\n"; +print $RCMD "plot_step_table(mygwas, \"extBIC\")\n"; +print $RCMD "plot_step_table(mygwas, \"maxpval\")\n"; +print $RCMD "plot_step_RSS(mygwas)\n"; +# for (my $j = 1; $j <= ($max_steps - 1); $j++) +# { + # print $RCMD "plot_fwd_GWAS(mygwas, step = $j, snp_info = map, pval_filt = 0.1)\n"; +# } +print $RCMD "plot_opt_GWAS(mygwas, opt = \"extBIC\", snp_info = map, pval_filt = 0.1)\n"; +print $RCMD "plot_opt_GWAS(mygwas, opt = \"mbonf\", snp_info = map, pval_filt = 0.1)\n"; +#print $RCMD "qqplot_fwd_GWAS(mygwas, nsteps = mlmm_data\$maxsteps)\n"; +print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"extBIC\")\n"; +print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"mbonf\")\n"; + +# outputs +print $RCMD "write.table(mygwas\$RSSout, '$out.rss', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n"; +print $RCMD "write.table(mygwas\$step_table, '$out.steptable', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n"; + +$plot_opt = "\$opt_" . $plot_opt; +print $RCMD "pval = mygwas" . $plot_opt . "\$out\n"; +print $RCMD "colnames(pval) = c(\"Marker_name\", \"Pvalue\")\n"; +print $RCMD "info_tmp = map\n"; +print $RCMD "colnames(info_tmp) = c(\"Marker_name\", \"Chr\", \"Pos\")\n"; +print $RCMD "res_asso = pval\n"; +print $RCMD qq~ +if(exists("info_tmp")){ + res_asso = merge(info_tmp, res_asso, by="Marker_name") + if( !is.element("Trait", colnames(info_tmp)) ){ + m = matrix(data="traitname", ncol=1, nrow=nrow(res_asso), dimnames=list(c(), c("Trait"))) + res_asso = cbind(m, res_asso) + } +} +~; +print $RCMD "res_asso = res_asso[order(res_asso[, \"Trait\"], res_asso[, \"Chr\"], res_asso[, \"Pos\"]), ]\n"; +print $RCMD "write.table(res_asso, '$out.res_asso', sep='\t', dec='.', quote=F, col.names=T, row.names=F)\n"; +close($RCMD); + +system("$RSCRIPT_EXE --vanilla rscript"); + + + + |
| b |
| diff -r 6b7107812931 -r 380b364980f9 MLMM.sh --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MLMM.sh Mon Apr 16 08:50:05 2018 -0400 |
| b |
| @@ -0,0 +1,28 @@ +#!/bin/bash +geno=$1 +map=$2 +pheno=$3 +steps=$4 +method=$5 +output=$6 +pdf=$7 +kinship=$8 +rss=$9 +step_table=${10} +log=${11} + +directory=`dirname $0` +mkdir tmpdir$$ +cp -rf $geno tmpdir$$/geno +cp -rf $map tmpdir$$/map +cp -rf $pheno tmpdir$$/pheno + +perl $directory/MLMM.pl -g tmpdir$$/geno -i tmpdir$$/map -p tmpdir$$/pheno -s $steps -m $method -o tmpdir$$/output -d $directory/source_library >>$log 2>&1 + +mv tmpdir$$/output.pdf $pdf +mv tmpdir$$/output.kinship $kinship +mv tmpdir$$/output.res_asso $output +mv tmpdir$$/output.rss $rss +mv tmpdir$$/output.steptable $step_table + + |
| b |
| diff -r 6b7107812931 -r 380b364980f9 MLMM.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MLMM.xml Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| @@ -0,0 +1,89 @@ +<tool id="mlmm" name="MLMM (GWAS analysis)" version="2.0.0"> + <description>GWAS using Multi-Locus Mixed-Model (MLMM)</description> + <requirements> + <requirement type="binary">Rscript</requirement> + <requirement type="binary">perl</requirement> + <requirement type="package" version="1.6.924">perl-bioperl</requirement> + </requirements> + <stdio> + <exit_code range="1:" level="fatal" /> + </stdio> + <command interpreter="bash">./MLMM.sh $geno $map $pheno $steps $method $output $pdf $kinship $rss $step_table $log</command> + <inputs> + <param format="txt" name="geno" type="data" label="Genotype matrix" help="NxM, N = individuals in line, M = Markers in columns, Genotype coded in 0,1,2"/> + <param type="data" format="txt" name="map" label="SNP Information file" help="3 columns: SNP, Chrom, Pos"/> + <param format="txt" name="pheno" type="data" label="Phenotype matrix" help="NxT, N = individuals in line, T = Trait in columns (Phenot1, Phenot2...)"/> + <param type="text" name="steps" label="Maximum number of steps for the forward approach" value="10"/> + <param name="method" type="select"> + <option value="extBIC">EBIC</option> + <option value="mbonf" selected="True">MBonf</option> + </param> + </inputs> + <outputs> + <data format="txt" name="output" label="Association results"/> + <data format="txt" name="kinship" label="Kinship matrix"/> + <data format="pdf" name="pdf" label="PDF Graphical outputs"/> + <data format="txt" name="rss" label="RSS"/> + <data format="txt" name="step_table" label="Step Table"/> + <data format="txt" name="log" label="Log file"/> + </outputs> + <tests> + <test> + <param name="geno" value="genot.txt"/> + <param name="map" value="map.txt"/> + <param name="pheno" value="phenot.txt"/> + <param name="steps" value="10"/> + <param name="method" value="mbonf"/> + <output name="output" value="output.txt"/> + <output name="kinship" value="kinship.txt"/> + <output name="rss" value="rss.txt"/> + <output name="step_table" value="step_table.txt"/> + </test> + </tests> + <help><![CDATA[ + + +.. class:: infomark + +**MLMM** version 1.0 + | **Please cite** "An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations.", **Segura V, Vilhjalmsson BJ, Platt A, Korte A, Seren U, Long Q, Nordborg M.**, Nature Genetics, 44: 825-830, 2012. + +.. class:: infomark + +**Galaxy integration** Provided by Southgreen & Dereeper Alexis (IRD) & Marcon Valentin (IFB & INRA) + +.. class:: infomark + +**Support** For any questions about Galaxy integration, please send an e-mail to alexis.dereeper@ird.fr + +--------------------------------------------------- + +==== +MLMM +==== + +----------- +Description +----------- + + | MLMM is an efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. + | For further informations, please visite the MLMM website_. + + +.. _website: https://sites.google.com/site/vincentosegura/mlmm + +------------ +Dependencies +------------ +MLMM + mlmm 1.0, Based on a archive provide by Vincent Segura: https://sites.google.com/site/vincentosegura/mlmm +Bioperl + perl-bioperl_ 1.6.924, Conda version + +.. _perl-bioperl: https://anaconda.org/bioconda/perl-bioperl + + ]]></help> + <citations> + <citation type="doi">10.1038/ng.2314</citation> + </citations> +</tool> |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/MLMM.pl --- a/mlmm/MLMM.pl Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,155 +0,0 @@ -#!/usr/bin/perl - -use strict; -use Switch; -use Getopt::Long; -use Bio::SeqIO; - - -my $usage = qq~Usage:$0 <args> [<opts>] -where <args> are: - -g, --geno <Genotype input> - -i, --info <SNP information. Genome position.> - -p, --pheno <Phenotype input> - -o, --out <output name> - -d, --directory <directory for MLMM R libraries> - -s, --step_number <number of steps. Maximum: 20. Default: 10> - -m, --method <Method: mbonf or extBIC. Default: mbonf> -~; -$usage .= "\n"; - -my ($geno,$map,$pheno,$out,$dir,$steps,$method); - - -GetOptions( - "geno=s" => \$geno, - "info=s" => \$map, - "pheno=s" => \$pheno, - "out=s" => \$out, - "dir=s" => \$dir, - "steps=s" => \$steps, - "method=s" => \$method -); - - -die $usage - if ( !$geno || !$map || !$pheno || !$out || !$dir || !$steps || !$method); - -my $max_steps = 10; -my $plot_opt = "mbonf"; -if ($method && $method ne 'mbonf' && $method ne 'extBIC') -{ - print "Aborted: Method must be mbonf or extBIC.\n"; - exit; -} -else -{ - $plot_opt = $method; -} -if ($steps && $steps !~/\d+/ && $steps > 20 && $steps < 2) -{ - print "Aborted: Number of steps must be greater than 2 and lower than 20.\n"; - exit; -} -else -{ - $max_steps = $steps; -} - - -my $chunk = 2; - -my $RSCRIPT_EXE = "Rscript"; -my $R_DIR = $dir; - - -my $head_trait = `head -1 $pheno`; -my @headers_traits = split(/\t/,$head_trait); -my $trait_name = $headers_traits[1]; - - -open( my $RCMD, ">rscript" ) or throw Error::Simple($!); - - - - -print $RCMD "Y_file <- \"" . $pheno . "\"\n"; -print $RCMD "X_file <- \"" . $geno . "\"\n"; -if($map) -{ - print $RCMD "map_file <- \"$map\"\n"; - print $RCMD "map <- read.table(map_file, sep = \"\\t\", header = T)\n"; -} -print $RCMD "mlmm_data = list()\n"; -print $RCMD "mlmm_data\$chunk <- $chunk\n"; -print $RCMD "mlmm_data\$maxsteps <- $max_steps\n"; -print $RCMD "genot <- read.table(X_file, sep = \"\\t\", header = T)\n"; -print $RCMD "genot_mat <- as.matrix(genot[, 2:ncol(genot)])\n"; -print $RCMD "rownames(genot_mat) <- genot\$Ind_id\n"; - -print $RCMD "phenot <- read.table(Y_file, sep = \"\\t\", header = T)\n"; - - - -# missing data imputation -print $RCMD "genot_imp <- genot_mat\n"; -print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n"; -print $RCMD "for (i in 1:ncol(genot_imp)){genot_imp[is.na(genot_imp[,i]), i] <- average[i]}\n"; - -# kinship matrix computation -print $RCMD "average <- colMeans(genot_imp, na.rm = T)\n"; -print $RCMD "stdev <- apply(genot_imp, 2, sd)\n"; -print $RCMD "genot_stand <- sweep(sweep(genot_imp, 2, average, \"-\"), 2, stdev, \"/\")\n"; -print $RCMD "K_mat <- (genot_stand %*% t(genot_stand)) / ncol(genot_stand)\n"; -print $RCMD "write.table(K_mat, '$out.kinship', sep='\\t', dec='.', quote=F, col.names=T, row.names=T)\n"; - -print $RCMD "source(\"" . $R_DIR. "/mlmm.r\")\n"; -print $RCMD "source(\"" . $R_DIR. "/emma.r\")\n"; - -# mlmm -print $RCMD "mygwas <- mlmm(Y = phenot\$$trait_name, X = genot_imp, K = K_mat, nbchunks=mlmm_data\$chunk, maxsteps=mlmm_data\$maxsteps)\n"; - -# plots -print $RCMD "pdf('$out.pdf')\n"; -print $RCMD "plot_step_table(mygwas, \"h2\")\n"; -print $RCMD "plot_step_table(mygwas, \"extBIC\")\n"; -print $RCMD "plot_step_table(mygwas, \"maxpval\")\n"; -print $RCMD "plot_step_RSS(mygwas)\n"; -# for (my $j = 1; $j <= ($max_steps - 1); $j++) -# { - # print $RCMD "plot_fwd_GWAS(mygwas, step = $j, snp_info = map, pval_filt = 0.1)\n"; -# } -print $RCMD "plot_opt_GWAS(mygwas, opt = \"extBIC\", snp_info = map, pval_filt = 0.1)\n"; -print $RCMD "plot_opt_GWAS(mygwas, opt = \"mbonf\", snp_info = map, pval_filt = 0.1)\n"; -#print $RCMD "qqplot_fwd_GWAS(mygwas, nsteps = mlmm_data\$maxsteps)\n"; -print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"extBIC\")\n"; -print $RCMD "qqplot_opt_GWAS(mygwas, opt = \"mbonf\")\n"; - -# outputs -print $RCMD "write.table(mygwas\$RSSout, '$out.rss', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n"; -print $RCMD "write.table(mygwas\$step_table, '$out.steptable', sep='\\t', dec='.', quote=F, col.names=T, row.names=F)\n"; - -$plot_opt = "\$opt_" . $plot_opt; -print $RCMD "pval = mygwas" . $plot_opt . "\$out\n"; -print $RCMD "colnames(pval) = c(\"Marker_name\", \"Pvalue\")\n"; -print $RCMD "info_tmp = map\n"; -print $RCMD "colnames(info_tmp) = c(\"Marker_name\", \"Chr\", \"Pos\")\n"; -print $RCMD "res_asso = pval\n"; -print $RCMD qq~ -if(exists("info_tmp")){ - res_asso = merge(info_tmp, res_asso, by="Marker_name") - if( !is.element("Trait", colnames(info_tmp)) ){ - m = matrix(data="traitname", ncol=1, nrow=nrow(res_asso), dimnames=list(c(), c("Trait"))) - res_asso = cbind(m, res_asso) - } -} -~; -print $RCMD "res_asso = res_asso[order(res_asso[, \"Trait\"], res_asso[, \"Chr\"], res_asso[, \"Pos\"]), ]\n"; -print $RCMD "write.table(res_asso, '$out.res_asso', sep='\t', dec='.', quote=F, col.names=T, row.names=F)\n"; -close($RCMD); - -system("$RSCRIPT_EXE --vanilla rscript"); - - - - |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/MLMM.sh --- a/mlmm/MLMM.sh Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,28 +0,0 @@ -#!/bin/bash -geno=$1 -map=$2 -pheno=$3 -steps=$4 -method=$5 -output=$6 -pdf=$7 -kinship=$8 -rss=$9 -step_table=${10} -log=${11} - -directory=`dirname $0` -mkdir tmpdir$$ -cp -rf $geno tmpdir$$/geno -cp -rf $map tmpdir$$/map -cp -rf $pheno tmpdir$$/pheno - -perl $directory/MLMM.pl -g tmpdir$$/geno -i tmpdir$$/map -p tmpdir$$/pheno -s $steps -m $method -o tmpdir$$/output -d $directory/source_library >>$log 2>&1 - -mv tmpdir$$/output.pdf $pdf -mv tmpdir$$/output.kinship $kinship -mv tmpdir$$/output.res_asso $output -mv tmpdir$$/output.rss $rss -mv tmpdir$$/output.steptable $step_table - - |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/MLMM.xml --- a/mlmm/MLMM.xml Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,61 +0,0 @@ -<tool id="mlmm" name="MLMM" version="1.0"> - <description>GWAS using Multi-Locus Mixed-Model (MLMM)</description> - <requirements> - <requirement type="binary">Rscript</requirement> - </requirements> - <command interpreter="bash">./MLMM.sh $geno $map $pheno $steps $method $output $pdf $kinship $rss $step_table $log - </command> - <inputs> - <param format="txt" name="geno" type="data" label="Genotype matrix" help="NxM, N = individuals in line, M = Markers in columns, Genotype coded in 0,1,2"/> - <param type="data" format="txt" name="map" label="SNP Information file" help="3 columns: SNP, Chrom, Pos"/> - <param format="txt" name="pheno" type="data" label="Phenotype matrix" help="NxT, N = individuals in line, T = Trait in columns (Phenot1, Phenot2...)"/> - <param type="text" name="steps" label="Maximum number of steps for the forward approach" value="10"/> - <param name="method" type="select"> - <option value="extBIC">EBIC</option> - <option value="mbonf" selected="True">MBonf</option> - </param> - </inputs> - <outputs> - <data format="txt" name="output" label="Association results"/> - <data format="txt" name="kinship" label="Kinship matrix"/> - <data format="pdf" name="pdf" label="PDF Graphical outputs"/> - <data format="txt" name="rss" label="RSS"/> - <data format="txt" name="step_table" label="Step Table"/> - <data format="txt" name="log" label="Log file"/> - </outputs> - <help> - - -.. class:: infomark - -**Program encapsulated in Galaxy by Southgreen** - -.. class:: infomark - -**MLMM version 1.0** - ------ - -============== - Please cite: -============== - -"An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations.", **Segura V, Vilhjalmsson BJ, Platt A, Korte A, Seren U, Long Q, Nordborg M.**, Nature Genetics, 44: 825-830, 2012. - ------ - -=========== - Overview: -=========== - -MLMM is an efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. - ------ - -For further informations, please visite the MLMM_ website. - - -.. _MLMM: https://sites.google.com/site/vincentosegura/mlmm - </help> - -</tool> |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/source_library/emma.mlmm.r --- a/mlmm/source_library/emma.mlmm.r Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,1274 +0,0 @@\n-emma.kinship <- function(snps, method="additive", use="all") {\n- n0 <- sum(snps==0,na.rm=TRUE)\n- nh <- sum(snps==0.5,na.rm=TRUE)\n- n1 <- sum(snps==1,na.rm=TRUE)\n- nNA <- sum(is.na(snps))\n-\n- stopifnot(n0+nh+n1+nNA == length(snps))\n-\n- if ( method == "dominant" ) {\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n- snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n- }\n- else if ( method == "recessive" ) {\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n- snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n- }\n- else if ( ( method == "additive" ) && ( nh > 0 ) ) {\n- dsnps <- snps\n- rsnps <- snps\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n- dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n- rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n- snps <- rbind(dsnps,rsnps)\n- }\n-\n- if ( use == "all" ) {\n- mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps))\n- snps[is.na(snps)] <- mafs[is.na(snps)]\n- }\n- else if ( use == "complete.obs" ) {\n- snps <- snps[rowSums(is.na(snps))==0,]\n- }\n-\n- n <- ncol(snps)\n- K <- matrix(nrow=n,ncol=n)\n- diag(K) <- 1\n-\n- for(i in 2:n) {\n- for(j in 1:(i-1)) {\n- x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j])\n- K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x))\n- K[j,i] <- K[i,j]\n- }\n- }\n- return(K)\n-}\n-\n-emma.eigen.L <- function(Z,K,complete=TRUE) {\n- if ( is.null(Z) ) {\n- return(emma.eigen.L.wo.Z(K))\n- }\n- else {\n- return(emma.eigen.L.w.Z(Z,K,complete))\n- }\n-}\n-\n-emma.eigen.L.wo.Z <- function(K) {\n- eig <- eigen(K,symmetric=TRUE)\n- return(list(values=eig$values,vectors=eig$vectors))\n-}\n-\n-emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {\n- if ( complete == FALSE ) {\n- vids <- colSums(Z)>0\n- Z <- Z[,vids]\n- K <- K[vids,vids]\n- }\n- eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)\n- return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))\n-}\n-\n-emma.eigen.R <- function(Z,K,X,complete=TRUE) {\n- if ( ncol(X) == 0 ) {\n- return(emma.eigen.L(Z,K))\n- }\n- else if ( is.null(Z) ) {\n- return(emma.eigen.R.wo.Z(K,X))\n- }\n- else {\n- return(emma.eigen.R.w.Z(Z,K,X,complete))\n- }\n-}\n-\n-emma.eigen.R.wo.Z <- function(K, X) {\n- n <- nrow(X)\n- q <- ncol(X)\n- S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)\n- eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)\n- stopifnot(!is.complex(eig$values))\n- return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))\n-}\n-\n-emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {\n- if ( complete == FALSE ) {\n- vids <- colSums(Z) > 0\n- Z <- Z[,vids]\n- K <- K[vids,vids]\n- }\n- n <- nrow(Z)\n- t <- ncol(Z)\n- q <- ncol(X)\n- \n- SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)\n- eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE)\n- if ( is.complex(eig$values) ) {\n- eig$values <- Re(eig$values)\n- eig$vectors <- Re(eig$vectors) \n- }\n- qr.X <- qr.Q(qr(X))\n- return(list(values=eig$values[1:(t-q)],\n- vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),\n- complete=TRUE)[,c(1:(t-q),(t+1):n)])) \n-}\n-\n-emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {\n- n <- length(xi)\n- delta <- exp(logdelta)\n- return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) ) \n-}\n-\n-emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {\n- t <- length(xi.1)\n- delta <- exp(logdelta)\n-# stopifnot(length(lambda) == length(etas.1))\n- return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) )\n-}\n-\n-emma.delta.ML.dLL.wo.Z <- functio'..b',drop=FALSE]))\n- if ( nv == t ) {\n- eig.R1 = emma.eigen.R.w.Z(Z,K,X)\n- } \n- }\n-\n- for(j in 1:g) {\n- vrows <- !is.na(ys[j,])\n- if ( nv == t ) {\n- yv <- ys[j,vrows]\n- nr <- sum(vrows)\n- if ( is.null(Z) ) {\n- if ( nr == n ) {\n- REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)\n- U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) \n- }\n- else {\n- eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])\n- REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n- }\n- dfs[i,j] <- nr-q1\n- }\n- else {\n- if ( nr == n ) {\n- REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)\n- U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) \n- }\n- else {\n- vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))\n- eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n- REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n- }\n- dfs[i,j] <- nr-q1\n- }\n-\n- yt <- crossprod(U,yv)\n- Xt <- crossprod(U,X[vrows,,drop=FALSE])\n- iXX <- solve(crossprod(Xt,Xt))\n- beta <- iXX%*%crossprod(Xt,yt)\n- if ( !ponly ) {\n- vgs[i,j] <- REMLE$vg\n- ves[i,j] <- REMLE$ve\n- REMLs[i,j] <- REMLE$REML\n- }\n- stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n- }\n- else {\n- if ( is.null(Z) ) {\n- vtids <- vrows & vids\n- eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])\n- yv <- ys[j,vtids]\n- nr <- sum(vtids)\n- REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n- Xt <- crossprod(U,X[vtids,,drop=FALSE])\n- dfs[i,j] <- nr-q1\n- }\n- else {\n- vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids\n- vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))\n- eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n- yv <- ys[j,vtrows]\n- nr <- sum(vtrows)\n- REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n- Xt <- crossprod(U,X[vtrows,,drop=FALSE])\n- dfs[i,j] <- nr-q1\n- }\n- yt <- crossprod(U,yv)\n- iXX <- solve(crossprod(Xt,Xt))\n- beta <- iXX%*%crossprod(Xt,yt)\n- if ( !ponly ) {\n- vgs[i,j] <- REMLE$vg\n- ves[i,j] <- REMLE$ve\n- REMLs[i,j] <- REMLE$REML\n- }\n- stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n- \n- }\n- }\n- ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) \n- }\n- } \n- }\n- if ( ponly ) {\n- return (ps)\n- }\n- else {\n- return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))\n- }\n-}\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/source_library/emma.r --- a/mlmm/source_library/emma.r Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,1274 +0,0 @@\n-emma.kinship <- function(snps, method="additive", use="all") {\n- n0 <- sum(snps==0,na.rm=TRUE)\n- nh <- sum(snps==0.5,na.rm=TRUE)\n- n1 <- sum(snps==1,na.rm=TRUE)\n- nNA <- sum(is.na(snps))\n-\n- stopifnot(n0+nh+n1+nNA == length(snps))\n-\n- if ( method == "dominant" ) {\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n- snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n- }\n- else if ( method == "recessive" ) {\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n- snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n- }\n- else if ( ( method == "additive" ) && ( nh > 0 ) ) {\n- dsnps <- snps\n- rsnps <- snps\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n- dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n- flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n- rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n- snps <- rbind(dsnps,rsnps)\n- }\n-\n- if ( use == "all" ) {\n- mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps))\n- snps[is.na(snps)] <- mafs[is.na(snps)]\n- }\n- else if ( use == "complete.obs" ) {\n- snps <- snps[rowSums(is.na(snps))==0,]\n- }\n-\n- n <- ncol(snps)\n- K <- matrix(nrow=n,ncol=n)\n- diag(K) <- 1\n-\n- for(i in 2:n) {\n- for(j in 1:(i-1)) {\n- x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j])\n- K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x))\n- K[j,i] <- K[i,j]\n- }\n- }\n- return(K)\n-}\n-\n-emma.eigen.L <- function(Z,K,complete=TRUE) {\n- if ( is.null(Z) ) {\n- return(emma.eigen.L.wo.Z(K))\n- }\n- else {\n- return(emma.eigen.L.w.Z(Z,K,complete))\n- }\n-}\n-\n-emma.eigen.L.wo.Z <- function(K) {\n- eig <- eigen(K,symmetric=TRUE)\n- return(list(values=eig$values,vectors=eig$vectors))\n-}\n-\n-emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {\n- if ( complete == FALSE ) {\n- vids <- colSums(Z)>0\n- Z <- Z[,vids]\n- K <- K[vids,vids]\n- }\n- eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)\n- return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))\n-}\n-\n-emma.eigen.R <- function(Z,K,X,complete=TRUE) {\n- if ( ncol(X) == 0 ) {\n- return(emma.eigen.L(Z,K))\n- }\n- else if ( is.null(Z) ) {\n- return(emma.eigen.R.wo.Z(K,X))\n- }\n- else {\n- return(emma.eigen.R.w.Z(Z,K,X,complete))\n- }\n-}\n-\n-emma.eigen.R.wo.Z <- function(K, X) {\n- n <- nrow(X)\n- q <- ncol(X)\n- S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)\n- eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)\n- stopifnot(!is.complex(eig$values))\n- return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))\n-}\n-\n-emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {\n- if ( complete == FALSE ) {\n- vids <- colSums(Z) > 0\n- Z <- Z[,vids]\n- K <- K[vids,vids]\n- }\n- n <- nrow(Z)\n- t <- ncol(Z)\n- q <- ncol(X)\n- \n- SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)\n- eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE)\n- if ( is.complex(eig$values) ) {\n- eig$values <- Re(eig$values)\n- eig$vectors <- Re(eig$vectors) \n- }\n- qr.X <- qr.Q(qr(X))\n- return(list(values=eig$values[1:(t-q)],\n- vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),\n- complete=TRUE)[,c(1:(t-q),(t+1):n)])) \n-}\n-\n-emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {\n- n <- length(xi)\n- delta <- exp(logdelta)\n- return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) ) \n-}\n-\n-emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {\n- t <- length(xi.1)\n- delta <- exp(logdelta)\n-# stopifnot(length(lambda) == length(etas.1))\n- return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) )\n-}\n-\n-emma.delta.ML.dLL.wo.Z <- functio'..b',drop=FALSE]))\n- if ( nv == t ) {\n- eig.R1 = emma.eigen.R.w.Z(Z,K,X)\n- } \n- }\n-\n- for(j in 1:g) {\n- vrows <- !is.na(ys[j,])\n- if ( nv == t ) {\n- yv <- ys[j,vrows]\n- nr <- sum(vrows)\n- if ( is.null(Z) ) {\n- if ( nr == n ) {\n- REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)\n- U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) \n- }\n- else {\n- eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])\n- REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n- }\n- dfs[i,j] <- nr-q1\n- }\n- else {\n- if ( nr == n ) {\n- REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)\n- U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) \n- }\n- else {\n- vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))\n- eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n- REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n- }\n- dfs[i,j] <- nr-q1\n- }\n-\n- yt <- crossprod(U,yv)\n- Xt <- crossprod(U,X[vrows,,drop=FALSE])\n- iXX <- solve(crossprod(Xt,Xt))\n- beta <- iXX%*%crossprod(Xt,yt)\n- if ( !ponly ) {\n- vgs[i,j] <- REMLE$vg\n- ves[i,j] <- REMLE$ve\n- REMLs[i,j] <- REMLE$REML\n- }\n- stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n- }\n- else {\n- if ( is.null(Z) ) {\n- vtids <- vrows & vids\n- eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])\n- yv <- ys[j,vtids]\n- nr <- sum(vtids)\n- REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n- Xt <- crossprod(U,X[vtids,,drop=FALSE])\n- dfs[i,j] <- nr-q1\n- }\n- else {\n- vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids\n- vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))\n- eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n- yv <- ys[j,vtrows]\n- nr <- sum(vtrows)\n- REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n- U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n- Xt <- crossprod(U,X[vtrows,,drop=FALSE])\n- dfs[i,j] <- nr-q1\n- }\n- yt <- crossprod(U,yv)\n- iXX <- solve(crossprod(Xt,Xt))\n- beta <- iXX%*%crossprod(Xt,yt)\n- if ( !ponly ) {\n- vgs[i,j] <- REMLE$vg\n- ves[i,j] <- REMLE$ve\n- REMLs[i,j] <- REMLE$REML\n- }\n- stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n- \n- }\n- }\n- ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) \n- }\n- } \n- }\n- if ( ponly ) {\n- return (ps)\n- }\n- else {\n- return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))\n- }\n-}\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/source_library/mlmm.r --- a/mlmm/source_library/mlmm.r Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b"@@ -1,477 +0,0 @@\n-##############################################################################################################################################\r\n-###MLMM - Multi-Locus Mixed Model\r\n-###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH\r\n-#######\r\n-#\r\n-##note: require EMMA\r\n-#library(emma)\r\n-#source('emma.r')\r\n-#\r\n-##REQUIRED DATA & FORMAT\r\n-#\r\n-#PHENOTYPE - Y: a vector of length m, with names(Y)=individual names\r\n-#GENOTYPE - X: a n by m matrix, where n=number of individuals, m=number of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names\r\n-#KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names\r\n-#each of these data being sorted in the same way, according to the individual name\r\n-#\r\n-##FOR PLOTING THE GWAS RESULTS\r\n-#SNP INFORMATION - snp_info: a data frame having at least 3 columns:\r\n-# - 1 named 'SNP', with SNP names (same as colnames(X)),\r\n-# - 1 named 'Chr', with the chromosome number to which belong each SNP\r\n-# - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.\r\n-#######\r\n-#\r\n-##FUNCTIONS USE\r\n-#save this file somewhere on your computer and source it!\r\n-#source('path/mlmm.r')\r\n-#\r\n-###FORWARD + BACKWARD ANALYSES\r\n-#mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)\r\n-#X,Y,K as described above\r\n-#nbchunks: an integer defining the number of chunks of X to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory\r\n-#maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,\r\n-#\t\t\thowever to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.\r\n-#\t\t\tIt's value must be specified as an integer >= 3\r\n-#\r\n-###RESULTS\r\n-#\r\n-##STEPWISE TABLE\r\n-#mygwas$step_table\r\n-#\r\n-##PLOTS\r\n-#\r\n-##PLOTS FORM THE FORWARD TABLE\r\n-#plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))\r\n-#\r\n-##RSS PLOT\r\n-#plot_step_RSS(mygwas)\r\n-#\r\n-##GWAS MANHATTAN PLOTS\r\n-#\r\n-#FORWARD STEPS\r\n-#plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)\r\n-#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#\r\n-#OPTIMAL MODELS\r\n-#Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria\r\n-#\r\n-#plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#\r\n-##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST\r\n-#plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)\r\n-#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#chrom is an integer specifying the chromosome on which the region of interest is\r\n-#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n-#\r\n-#plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#chrom is an integer specifying the chromosome on which the region of interest is\r\n-#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n-#\r\n-##QQPLOTS of pvalues\r\n-#qqplot_fwd_GWAS(mygwas,nsteps)\r\n-#nsteps=maximum number of forward steps to be displayed\r\n-#\r\n-#qqplot_opt_GWAS(mygwas,opt=c('extBIC"..b'ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2\r\n-\r\n-\tplot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab=\'-log10(pval)\',xaxt=\'n\',xlab=\'chromosome\')\r\n-\taxis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))\r\n-\tabline(h=x$bonf_thresh,lty=3,col=\'black\')}\r\n-\r\n-plot_region<-function(x,chrom,pos1,pos2){\r\n-\tregion<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)\r\n-\tregion$col<- if (chrom %% 2 == 0) {\'gray60\'} else {\'gray10\'}\r\n-\tregion$col[which(region$SNP %in% x$cof)]<-\'red\'\r\n-\tplot(region$Pos,-log10(region$pval),type=\'p\',pch=20,main=paste(\'chromosome\',chrom,sep=\'\'),xlab=\'position (bp)\',ylab=\'-log10(pval)\',col=region$col,xlim=c(pos1,pos2))\r\n-\tabline(h=x$bonf_thresh,lty=3,col=\'black\')}\r\n-\r\n-\r\n-plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {\r\n-\tstopifnot(step<=length(x$pval_step))\r\n-\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n-\tplot_GWAS(output)}\r\n-\r\n-plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {\r\n-\tstopifnot(step<=length(x$pval_step))\r\n-\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n-\tplot_region(output,chrom,pos1,pos2)}\r\n-\r\n-\r\n-plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {\r\n-\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_GWAS(output)}\r\n-\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_GWAS(output)}\r\n-\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n-\r\n-plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {\r\n-\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_region(output,chrom,pos1,pos2)}\r\n-\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_region(output,chrom,pos1,pos2)}\r\n-\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n-\r\n-\r\n-qqplot_fwd_GWAS<-function(x,nsteps){\r\n-\tstopifnot(nsteps<=length(x$pval_step))\r\n-\te<--log10(ppoints(nrow(x$pval_step[[1]]$out)))\r\n-\tostep<-list()\r\n-\tostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))\r\n-\tfor (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}\r\n-\r\n-\tmaxp<-ceiling(max(unlist(ostep)))\r\n-\r\n-\tplot(e,ostep[[1]],type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))\r\n-\tabline(0,1,col="dark grey")\r\n-\r\n-\tfor (i in 2:nsteps) {\r\n-\tpar(new=T)\r\n-\tplot(e,ostep[[i]],type=\'b\',pch=20,cex=0.8,col=i,axes=\'F\',xlab=\'\',ylab=\'\',xlim=c(0,max(e)+1),ylim=c(0,maxp))}\r\n-\tlegend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),\'cof\',sep=\' \'))\r\n-}\r\n-\r\n-qqplot_opt_GWAS<-function(x,opt){\r\n-\tif (opt==\'extBIC\') {\r\n-\t\te<--log10(ppoints(nrow(x$opt_extBIC$out)))\r\n-\t\to<--log10(sort(x$opt_extBIC$out$pval))\r\n-\t\tmaxp<-ceiling(max(o))\r\n-\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to extBIC\'))\r\n-\t\tabline(0,1,col="dark grey")}\r\n-\telse if (opt==\'mbonf\') {\r\n-\t\te<--log10(ppoints(nrow(x$opt_mbonf$out)))\r\n-\t\to<--log10(sort(x$opt_mbonf$out$pval))\r\n-\t\tmaxp<-ceiling(max(o))\r\n-\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to mbonf\'))\r\n-\t\tabline(0,1,col="dark grey")}\r\n-\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n-\r\n-\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/source_library/mlmm1.r --- a/mlmm/source_library/mlmm1.r Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b"@@ -1,486 +0,0 @@\n-##############################################################################################################################################\r\n-###MLMM - Multi-Locus Mixed Model\r\n-###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH\r\n-#######\r\n-#\r\n-##note: require EMMA\r\n-#library(emma)\r\n-#source('emma.r')\r\n-#\r\n-##REQUIRED DATA & FORMAT\r\n-#\r\n-#PHENOTYPE - Y: a vector of length m, with names(Y)=individual names\r\n-#GENOTYPE - X: a n by m matrix, where n=number of individuals, m=n umber of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names\r\n-#KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names\r\n-#each of these data being sorted in the same way, according to the individual name\r\n-#\r\n-##FOR PLOTING THE GWAS RESULTS\r\n-#SNP INFORMATION - snp_info: a data frame having at least 3 columns:\r\n-# - 1 named 'SNP', with SNP names (same as colnames(X)),\r\n-# - 1 named 'Chr', with the chromosome number to which belong each SNP\r\n-# - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.\r\n-#######\r\n-#\r\n-##FUNCTIONS USE\r\n-#save this file somewhere on your computer and source it!\r\n-#source('path/mlmm.r')\r\n-#\r\n-###FORWARD + BACKWARD ANALYSES\r\n-#mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)\r\n-#X,Y,K as described above\r\n-#nbchunks: an integer defining the number of chunks of to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory\r\n-#maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,\r\n-#\t\t\thowever to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.\r\n-#\t\t\tIt's value must be specified as an integer >= 3\r\n-#\r\n-###RESULTS\r\n-#\r\n-##STEPWISE TABLE\r\n-#mygwas$step_table\r\n-#\r\n-##PLOTS\r\n-#\r\n-##PLOTS FORM THE FORWARD TABLE\r\n-#plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))\r\n-#\r\n-##RSS PLOT\r\n-#plot_step_RSS(mygwas)\r\n-#\r\n-##GWAS MANHATTAN PLOTS\r\n-#\r\n-#FORWARD STEPS\r\n-#plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)\r\n-#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#\r\n-#OPTIMAL MODELS\r\n-#Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria\r\n-#\r\n-#plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#\r\n-##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST\r\n-#plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)\r\n-#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#chrom is an integer specifying the chromosome on which the region of interest is\r\n-#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n-#\r\n-#plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)\r\n-#snp_info as described above\r\n-#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n-#chrom is an integer specifying the chromosome on which the region of interest is\r\n-#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n-#\r\n-##QQPLOTS of pvalues\r\n-#qqplot_fwd_GWAS(mygwas,nsteps)\r\n-#nsteps=maximum number of forward steps to be displayed\r\n-#\r\n-#qqplot_opt_GWAS(mygwas,opt=c('extBIC'"..b'output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7) \r\n- legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))\r\n- } else {\r\n- legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))\r\n- }\r\n-}\r\n- \r\n-plot_region<-function(x,chrom,pos1,pos2){\r\n-\tregion<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)\r\n-\tregion$col<- if (chrom %% 2 == 0) {\'gray60\'} else {\'gray10\'}\r\n-\tregion$col[which(region$SNP %in% x$cof)]<-\'red\'\r\n-\tplot(region$Pos,-log10(region$pval),type=\'p\',pch=20,main=paste(\'chromosome\',chrom,sep=\'\'),xlab=\'position (bp)\',ylab=\'-log10(pval)\',col=region$col,xlim=c(pos1,pos2))\r\n-\tabline(h=x$bonf_thresh,lty=3,col=\'black\')}\r\n-\r\n-\r\n-plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {\r\n-\tstopifnot(step<=length(x$pval_step))\r\n-\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n-\tplot_GWAS(output)}\r\n-\r\n-plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {\r\n-\tstopifnot(step<=length(x$pval_step))\r\n-\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n-\tplot_region(output,chrom,pos1,pos2)}\r\n-\r\n-\r\n-plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {\r\n-\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_GWAS(output)}\r\n-\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_GWAS(output)}\r\n-\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n-\r\n-plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {\r\n-\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_region(output,chrom,pos1,pos2)}\r\n-\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n-\t\tplot_region(output,chrom,pos1,pos2)}\r\n-\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n-\r\n-\r\n-qqplot_fwd_GWAS<-function(x,nsteps){\r\n-\tstopifnot(nsteps<=length(x$pval_step))\r\n-\te<--log10(ppoints(nrow(x$pval_step[[1]]$out)))\r\n-\tostep<-list()\r\n-\tostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))\r\n-\tfor (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}\r\n-\r\n-\tmaxp<-ceiling(max(unlist(ostep)))\r\n-\r\n-\tplot(e,ostep[[1]],type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))\r\n-\tabline(0,1,col="dark grey")\r\n-\r\n-\tfor (i in 2:nsteps) {\r\n-\tpar(new=T)\r\n-\tplot(e,ostep[[i]],type=\'b\',pch=20,cex=0.8,col=i,axes=\'F\',xlab=\'\',ylab=\'\',xlim=c(0,max(e)+1),ylim=c(0,maxp))}\r\n-\tlegend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),\'cof\',sep=\' \'))\r\n-}\r\n-\r\n-qqplot_opt_GWAS<-function(x,opt){\r\n-\tif (opt==\'extBIC\') {\r\n-\t\te<--log10(ppoints(nrow(x$opt_extBIC$out)))\r\n-\t\to<--log10(sort(x$opt_extBIC$out$pval))\r\n-\t\tmaxp<-ceiling(max(o))\r\n-\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to extBIC\'))\r\n-\t\tabline(0,1,col="dark grey")}\r\n-\telse if (opt==\'mbonf\') {\r\n-\t\te<--log10(ppoints(nrow(x$opt_mbonf$out)))\r\n-\t\to<--log10(sort(x$opt_mbonf$out$pval))\r\n-\t\tmaxp<-ceiling(max(o))\r\n-\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to mbonf\'))\r\n-\t\tabline(0,1,col="dark grey")}\r\n-\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n-\r\n-\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 mlmm/source_library/plot_MLMM.Ago.r --- a/mlmm/source_library/plot_MLMM.Ago.r Thu Jul 02 05:42:38 2015 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,36 +0,0 @@ -## fonction ########################################################################### -plot_MLMM.Ago<-function(x) { - output1 <-x$out[order(x$out$Pos),] - output_ok<-output1[order(output1$Chr),] - - maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200))) - plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2)) -# plot_col<-c('blue','darkgreen','red','cyan','purple') - size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x - a<-rep(maxpos[1],size[1]) - b<-rep(plot_col[1],size[1]) - for (i in 2:max(unique(output_ok$Chr))){ - a<-c(a,rep(maxpos[i],size[i])) - b<-c(b,rep(plot_col[i],size[i]))} - - output_ok$xpos<-output_ok$Pos+a - output_ok$col<-b - output_ok$col[output_ok$mk=='qtl']<-'cyan' - output_ok$col[output_ok$SNP %in% x$cof]<-'red' - - d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2 - - plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome') - axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr)))) - xline(output_ok$xpos[output_ok$mk=='qtl'], col='cyan', lwd=0.1) - abline(h=x$bonf_thresh,lty=3,col='black') - - - if (length(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]) > 0) { - text(output_ok$xpos[-log10(output_ok$pval) > x$bonf_thresh], -log10(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7) - legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" ")) - } else { - legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" ")) - } - -} |
| b |
| diff -r 6b7107812931 -r 380b364980f9 source_library/emma.mlmm.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/source_library/emma.mlmm.r Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| b'@@ -0,0 +1,1274 @@\n+emma.kinship <- function(snps, method="additive", use="all") {\n+ n0 <- sum(snps==0,na.rm=TRUE)\n+ nh <- sum(snps==0.5,na.rm=TRUE)\n+ n1 <- sum(snps==1,na.rm=TRUE)\n+ nNA <- sum(is.na(snps))\n+\n+ stopifnot(n0+nh+n1+nNA == length(snps))\n+\n+ if ( method == "dominant" ) {\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n+ snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n+ }\n+ else if ( method == "recessive" ) {\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n+ snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n+ }\n+ else if ( ( method == "additive" ) && ( nh > 0 ) ) {\n+ dsnps <- snps\n+ rsnps <- snps\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n+ dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n+ rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n+ snps <- rbind(dsnps,rsnps)\n+ }\n+\n+ if ( use == "all" ) {\n+ mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps))\n+ snps[is.na(snps)] <- mafs[is.na(snps)]\n+ }\n+ else if ( use == "complete.obs" ) {\n+ snps <- snps[rowSums(is.na(snps))==0,]\n+ }\n+\n+ n <- ncol(snps)\n+ K <- matrix(nrow=n,ncol=n)\n+ diag(K) <- 1\n+\n+ for(i in 2:n) {\n+ for(j in 1:(i-1)) {\n+ x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j])\n+ K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x))\n+ K[j,i] <- K[i,j]\n+ }\n+ }\n+ return(K)\n+}\n+\n+emma.eigen.L <- function(Z,K,complete=TRUE) {\n+ if ( is.null(Z) ) {\n+ return(emma.eigen.L.wo.Z(K))\n+ }\n+ else {\n+ return(emma.eigen.L.w.Z(Z,K,complete))\n+ }\n+}\n+\n+emma.eigen.L.wo.Z <- function(K) {\n+ eig <- eigen(K,symmetric=TRUE)\n+ return(list(values=eig$values,vectors=eig$vectors))\n+}\n+\n+emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {\n+ if ( complete == FALSE ) {\n+ vids <- colSums(Z)>0\n+ Z <- Z[,vids]\n+ K <- K[vids,vids]\n+ }\n+ eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)\n+ return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))\n+}\n+\n+emma.eigen.R <- function(Z,K,X,complete=TRUE) {\n+ if ( ncol(X) == 0 ) {\n+ return(emma.eigen.L(Z,K))\n+ }\n+ else if ( is.null(Z) ) {\n+ return(emma.eigen.R.wo.Z(K,X))\n+ }\n+ else {\n+ return(emma.eigen.R.w.Z(Z,K,X,complete))\n+ }\n+}\n+\n+emma.eigen.R.wo.Z <- function(K, X) {\n+ n <- nrow(X)\n+ q <- ncol(X)\n+ S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)\n+ eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)\n+ stopifnot(!is.complex(eig$values))\n+ return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))\n+}\n+\n+emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {\n+ if ( complete == FALSE ) {\n+ vids <- colSums(Z) > 0\n+ Z <- Z[,vids]\n+ K <- K[vids,vids]\n+ }\n+ n <- nrow(Z)\n+ t <- ncol(Z)\n+ q <- ncol(X)\n+ \n+ SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)\n+ eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE)\n+ if ( is.complex(eig$values) ) {\n+ eig$values <- Re(eig$values)\n+ eig$vectors <- Re(eig$vectors) \n+ }\n+ qr.X <- qr.Q(qr(X))\n+ return(list(values=eig$values[1:(t-q)],\n+ vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),\n+ complete=TRUE)[,c(1:(t-q),(t+1):n)])) \n+}\n+\n+emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {\n+ n <- length(xi)\n+ delta <- exp(logdelta)\n+ return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) ) \n+}\n+\n+emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {\n+ t <- length(xi.1)\n+ delta <- exp(logdelta)\n+# stopifnot(length(lambda) == length(etas.1))\n+ return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) )\n+}\n+\n+emma.delta.ML.dLL.wo.Z <- functio'..b',drop=FALSE]))\n+ if ( nv == t ) {\n+ eig.R1 = emma.eigen.R.w.Z(Z,K,X)\n+ } \n+ }\n+\n+ for(j in 1:g) {\n+ vrows <- !is.na(ys[j,])\n+ if ( nv == t ) {\n+ yv <- ys[j,vrows]\n+ nr <- sum(vrows)\n+ if ( is.null(Z) ) {\n+ if ( nr == n ) {\n+ REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)\n+ U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) \n+ }\n+ else {\n+ eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])\n+ REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n+ }\n+ dfs[i,j] <- nr-q1\n+ }\n+ else {\n+ if ( nr == n ) {\n+ REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)\n+ U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) \n+ }\n+ else {\n+ vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))\n+ eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n+ REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n+ }\n+ dfs[i,j] <- nr-q1\n+ }\n+\n+ yt <- crossprod(U,yv)\n+ Xt <- crossprod(U,X[vrows,,drop=FALSE])\n+ iXX <- solve(crossprod(Xt,Xt))\n+ beta <- iXX%*%crossprod(Xt,yt)\n+ if ( !ponly ) {\n+ vgs[i,j] <- REMLE$vg\n+ ves[i,j] <- REMLE$ve\n+ REMLs[i,j] <- REMLE$REML\n+ }\n+ stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n+ }\n+ else {\n+ if ( is.null(Z) ) {\n+ vtids <- vrows & vids\n+ eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])\n+ yv <- ys[j,vtids]\n+ nr <- sum(vtids)\n+ REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n+ Xt <- crossprod(U,X[vtids,,drop=FALSE])\n+ dfs[i,j] <- nr-q1\n+ }\n+ else {\n+ vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids\n+ vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))\n+ eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n+ yv <- ys[j,vtrows]\n+ nr <- sum(vtrows)\n+ REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n+ Xt <- crossprod(U,X[vtrows,,drop=FALSE])\n+ dfs[i,j] <- nr-q1\n+ }\n+ yt <- crossprod(U,yv)\n+ iXX <- solve(crossprod(Xt,Xt))\n+ beta <- iXX%*%crossprod(Xt,yt)\n+ if ( !ponly ) {\n+ vgs[i,j] <- REMLE$vg\n+ ves[i,j] <- REMLE$ve\n+ REMLs[i,j] <- REMLE$REML\n+ }\n+ stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n+ \n+ }\n+ }\n+ ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) \n+ }\n+ } \n+ }\n+ if ( ponly ) {\n+ return (ps)\n+ }\n+ else {\n+ return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))\n+ }\n+}\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 source_library/emma.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/source_library/emma.r Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| b'@@ -0,0 +1,1274 @@\n+emma.kinship <- function(snps, method="additive", use="all") {\n+ n0 <- sum(snps==0,na.rm=TRUE)\n+ nh <- sum(snps==0.5,na.rm=TRUE)\n+ n1 <- sum(snps==1,na.rm=TRUE)\n+ nNA <- sum(is.na(snps))\n+\n+ stopifnot(n0+nh+n1+nNA == length(snps))\n+\n+ if ( method == "dominant" ) {\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n+ snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n+ }\n+ else if ( method == "recessive" ) {\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n+ snps[!is.na(snps) & (snps == 0.5)] <- flags[!is.na(snps) & (snps == 0.5)]\n+ }\n+ else if ( ( method == "additive" ) && ( nh > 0 ) ) {\n+ dsnps <- snps\n+ rsnps <- snps\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) > 0.5),nrow(snps),ncol(snps))\n+ dsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n+ flags <- matrix(as.double(rowMeans(snps,na.rm=TRUE) < 0.5),nrow(snps),ncol(snps))\n+ rsnps[!is.na(snps) & (snps==0.5)] <- flags[!is.na(snps) & (snps==0.5)]\n+ snps <- rbind(dsnps,rsnps)\n+ }\n+\n+ if ( use == "all" ) {\n+ mafs <- matrix(rowMeans(snps,na.rm=TRUE),nrow(snps),ncol(snps))\n+ snps[is.na(snps)] <- mafs[is.na(snps)]\n+ }\n+ else if ( use == "complete.obs" ) {\n+ snps <- snps[rowSums(is.na(snps))==0,]\n+ }\n+\n+ n <- ncol(snps)\n+ K <- matrix(nrow=n,ncol=n)\n+ diag(K) <- 1\n+\n+ for(i in 2:n) {\n+ for(j in 1:(i-1)) {\n+ x <- snps[,i]*snps[,j] + (1-snps[,i])*(1-snps[,j])\n+ K[i,j] <- sum(x,na.rm=TRUE)/sum(!is.na(x))\n+ K[j,i] <- K[i,j]\n+ }\n+ }\n+ return(K)\n+}\n+\n+emma.eigen.L <- function(Z,K,complete=TRUE) {\n+ if ( is.null(Z) ) {\n+ return(emma.eigen.L.wo.Z(K))\n+ }\n+ else {\n+ return(emma.eigen.L.w.Z(Z,K,complete))\n+ }\n+}\n+\n+emma.eigen.L.wo.Z <- function(K) {\n+ eig <- eigen(K,symmetric=TRUE)\n+ return(list(values=eig$values,vectors=eig$vectors))\n+}\n+\n+emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {\n+ if ( complete == FALSE ) {\n+ vids <- colSums(Z)>0\n+ Z <- Z[,vids]\n+ K <- K[vids,vids]\n+ }\n+ eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)\n+ return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))\n+}\n+\n+emma.eigen.R <- function(Z,K,X,complete=TRUE) {\n+ if ( ncol(X) == 0 ) {\n+ return(emma.eigen.L(Z,K))\n+ }\n+ else if ( is.null(Z) ) {\n+ return(emma.eigen.R.wo.Z(K,X))\n+ }\n+ else {\n+ return(emma.eigen.R.w.Z(Z,K,X,complete))\n+ }\n+}\n+\n+emma.eigen.R.wo.Z <- function(K, X) {\n+ n <- nrow(X)\n+ q <- ncol(X)\n+ S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)\n+ eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)\n+ stopifnot(!is.complex(eig$values))\n+ return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))\n+}\n+\n+emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {\n+ if ( complete == FALSE ) {\n+ vids <- colSums(Z) > 0\n+ Z <- Z[,vids]\n+ K <- K[vids,vids]\n+ }\n+ n <- nrow(Z)\n+ t <- ncol(Z)\n+ q <- ncol(X)\n+ \n+ SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)\n+ eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE,EISPACK=TRUE)\n+ if ( is.complex(eig$values) ) {\n+ eig$values <- Re(eig$values)\n+ eig$vectors <- Re(eig$vectors) \n+ }\n+ qr.X <- qr.Q(qr(X))\n+ return(list(values=eig$values[1:(t-q)],\n+ vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),\n+ complete=TRUE)[,c(1:(t-q),(t+1):n)])) \n+}\n+\n+emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {\n+ n <- length(xi)\n+ delta <- exp(logdelta)\n+ return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(lambda+delta))))-sum(log(xi+delta))) ) \n+}\n+\n+emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {\n+ t <- length(xi.1)\n+ delta <- exp(logdelta)\n+# stopifnot(length(lambda) == length(etas.1))\n+ return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(lambda+delta))+etas.2.sq/delta))-(sum(log(xi.1+delta))+(n-t)*logdelta)) )\n+}\n+\n+emma.delta.ML.dLL.wo.Z <- functio'..b',drop=FALSE]))\n+ if ( nv == t ) {\n+ eig.R1 = emma.eigen.R.w.Z(Z,K,X)\n+ } \n+ }\n+\n+ for(j in 1:g) {\n+ vrows <- !is.na(ys[j,])\n+ if ( nv == t ) {\n+ yv <- ys[j,vrows]\n+ nr <- sum(vrows)\n+ if ( is.null(Z) ) {\n+ if ( nr == n ) {\n+ REMLE <- emma.REMLE(yv,X,K,NULL,ngrids,llim,ulim,esp,eig.R1)\n+ U <- eig.L$vectors * matrix(sqrt(1/(eig.L$values+REMLE$delta)),n,n,byrow=TRUE) \n+ }\n+ else {\n+ eig.L0 <- emma.eigen.L.wo.Z(K[vrows,vrows,drop=FALSE])\n+ REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vrows,vrows,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n+ }\n+ dfs[i,j] <- nr-q1\n+ }\n+ else {\n+ if ( nr == n ) {\n+ REMLE <- emma.REMLE(yv,X,K,Z,ngrids,llim,ulim,esp,eig.R1)\n+ U <- eig.L$vectors * matrix(c(sqrt(1/(eig.L$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),n-t)),n,n,byrow=TRUE) \n+ }\n+ else {\n+ vtids <- as.logical(colSums(Z[vrows,,drop=FALSE]))\n+ eig.L0 <- emma.eigen.L.w.Z(Z[vrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n+ REMLE <- emma.REMLE(yv,X[vrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n+ }\n+ dfs[i,j] <- nr-q1\n+ }\n+\n+ yt <- crossprod(U,yv)\n+ Xt <- crossprod(U,X[vrows,,drop=FALSE])\n+ iXX <- solve(crossprod(Xt,Xt))\n+ beta <- iXX%*%crossprod(Xt,yt)\n+ if ( !ponly ) {\n+ vgs[i,j] <- REMLE$vg\n+ ves[i,j] <- REMLE$ve\n+ REMLs[i,j] <- REMLE$REML\n+ }\n+ stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n+ }\n+ else {\n+ if ( is.null(Z) ) {\n+ vtids <- vrows & vids\n+ eig.L0 <- emma.eigen.L.wo.Z(K[vtids,vtids,drop=FALSE])\n+ yv <- ys[j,vtids]\n+ nr <- sum(vtids)\n+ REMLE <- emma.REMLE(yv,X[vtids,,drop=FALSE],K[vtids,vtids,drop=FALSE],NULL,ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(sqrt(1/(eig.L0$values+REMLE$delta)),nr,nr,byrow=TRUE)\n+ Xt <- crossprod(U,X[vtids,,drop=FALSE])\n+ dfs[i,j] <- nr-q1\n+ }\n+ else {\n+ vtids <- as.logical(colSums(Z[vrows,,drop=FALSE])) & vids\n+ vtrows <- vrows & as.logical(rowSums(Z[,vids,drop=FALSE]))\n+ eig.L0 <- emma.eigen.L.w.Z(Z[vtrows,vtids,drop=FALSE],K[vtids,vtids,drop=FALSE])\n+ yv <- ys[j,vtrows]\n+ nr <- sum(vtrows)\n+ REMLE <- emma.REMLE(yv,X[vtrows,,drop=FALSE],K[vtids,vtids,drop=FALSE],Z[vtrows,vtids,drop=FALSE],ngrids,llim,ulim,esp)\n+ U <- eig.L0$vectors * matrix(c(sqrt(1/(eig.L0$values+REMLE$delta)),rep(sqrt(1/REMLE$delta),nr-sum(vtids))),nr,nr,byrow=TRUE)\n+ Xt <- crossprod(U,X[vtrows,,drop=FALSE])\n+ dfs[i,j] <- nr-q1\n+ }\n+ yt <- crossprod(U,yv)\n+ iXX <- solve(crossprod(Xt,Xt))\n+ beta <- iXX%*%crossprod(Xt,yt)\n+ if ( !ponly ) {\n+ vgs[i,j] <- REMLE$vg\n+ ves[i,j] <- REMLE$ve\n+ REMLs[i,j] <- REMLE$REML\n+ }\n+ stats[i,j] <- beta[q1]/sqrt(iXX[q1,q1]*REMLE$vg)\n+ \n+ }\n+ }\n+ ps[i,] <- 2*pt(abs(stats[i,]),dfs[i,],lower.tail=FALSE) \n+ }\n+ } \n+ }\n+ if ( ponly ) {\n+ return (ps)\n+ }\n+ else {\n+ return (list(ps=ps,REMLs=REMLs,stats=stats,dfs=dfs,vgs=vgs,ves=ves))\n+ }\n+}\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 source_library/mlmm.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/source_library/mlmm.r Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| b"@@ -0,0 +1,477 @@\n+##############################################################################################################################################\r\n+###MLMM - Multi-Locus Mixed Model\r\n+###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH\r\n+#######\r\n+#\r\n+##note: require EMMA\r\n+#library(emma)\r\n+#source('emma.r')\r\n+#\r\n+##REQUIRED DATA & FORMAT\r\n+#\r\n+#PHENOTYPE - Y: a vector of length m, with names(Y)=individual names\r\n+#GENOTYPE - X: a n by m matrix, where n=number of individuals, m=number of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names\r\n+#KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names\r\n+#each of these data being sorted in the same way, according to the individual name\r\n+#\r\n+##FOR PLOTING THE GWAS RESULTS\r\n+#SNP INFORMATION - snp_info: a data frame having at least 3 columns:\r\n+# - 1 named 'SNP', with SNP names (same as colnames(X)),\r\n+# - 1 named 'Chr', with the chromosome number to which belong each SNP\r\n+# - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.\r\n+#######\r\n+#\r\n+##FUNCTIONS USE\r\n+#save this file somewhere on your computer and source it!\r\n+#source('path/mlmm.r')\r\n+#\r\n+###FORWARD + BACKWARD ANALYSES\r\n+#mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)\r\n+#X,Y,K as described above\r\n+#nbchunks: an integer defining the number of chunks of X to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory\r\n+#maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,\r\n+#\t\t\thowever to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.\r\n+#\t\t\tIt's value must be specified as an integer >= 3\r\n+#\r\n+###RESULTS\r\n+#\r\n+##STEPWISE TABLE\r\n+#mygwas$step_table\r\n+#\r\n+##PLOTS\r\n+#\r\n+##PLOTS FORM THE FORWARD TABLE\r\n+#plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))\r\n+#\r\n+##RSS PLOT\r\n+#plot_step_RSS(mygwas)\r\n+#\r\n+##GWAS MANHATTAN PLOTS\r\n+#\r\n+#FORWARD STEPS\r\n+#plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)\r\n+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#\r\n+#OPTIMAL MODELS\r\n+#Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria\r\n+#\r\n+#plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#\r\n+##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST\r\n+#plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)\r\n+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#chrom is an integer specifying the chromosome on which the region of interest is\r\n+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n+#\r\n+#plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#chrom is an integer specifying the chromosome on which the region of interest is\r\n+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n+#\r\n+##QQPLOTS of pvalues\r\n+#qqplot_fwd_GWAS(mygwas,nsteps)\r\n+#nsteps=maximum number of forward steps to be displayed\r\n+#\r\n+#qqplot_opt_GWAS(mygwas,opt=c('extBIC"..b'ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2\r\n+\r\n+\tplot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab=\'-log10(pval)\',xaxt=\'n\',xlab=\'chromosome\')\r\n+\taxis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))\r\n+\tabline(h=x$bonf_thresh,lty=3,col=\'black\')}\r\n+\r\n+plot_region<-function(x,chrom,pos1,pos2){\r\n+\tregion<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)\r\n+\tregion$col<- if (chrom %% 2 == 0) {\'gray60\'} else {\'gray10\'}\r\n+\tregion$col[which(region$SNP %in% x$cof)]<-\'red\'\r\n+\tplot(region$Pos,-log10(region$pval),type=\'p\',pch=20,main=paste(\'chromosome\',chrom,sep=\'\'),xlab=\'position (bp)\',ylab=\'-log10(pval)\',col=region$col,xlim=c(pos1,pos2))\r\n+\tabline(h=x$bonf_thresh,lty=3,col=\'black\')}\r\n+\r\n+\r\n+plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {\r\n+\tstopifnot(step<=length(x$pval_step))\r\n+\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n+\tplot_GWAS(output)}\r\n+\r\n+plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {\r\n+\tstopifnot(step<=length(x$pval_step))\r\n+\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n+\tplot_region(output,chrom,pos1,pos2)}\r\n+\r\n+\r\n+plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {\r\n+\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_GWAS(output)}\r\n+\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_GWAS(output)}\r\n+\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n+\r\n+plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {\r\n+\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_region(output,chrom,pos1,pos2)}\r\n+\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_region(output,chrom,pos1,pos2)}\r\n+\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n+\r\n+\r\n+qqplot_fwd_GWAS<-function(x,nsteps){\r\n+\tstopifnot(nsteps<=length(x$pval_step))\r\n+\te<--log10(ppoints(nrow(x$pval_step[[1]]$out)))\r\n+\tostep<-list()\r\n+\tostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))\r\n+\tfor (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}\r\n+\r\n+\tmaxp<-ceiling(max(unlist(ostep)))\r\n+\r\n+\tplot(e,ostep[[1]],type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))\r\n+\tabline(0,1,col="dark grey")\r\n+\r\n+\tfor (i in 2:nsteps) {\r\n+\tpar(new=T)\r\n+\tplot(e,ostep[[i]],type=\'b\',pch=20,cex=0.8,col=i,axes=\'F\',xlab=\'\',ylab=\'\',xlim=c(0,max(e)+1),ylim=c(0,maxp))}\r\n+\tlegend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),\'cof\',sep=\' \'))\r\n+}\r\n+\r\n+qqplot_opt_GWAS<-function(x,opt){\r\n+\tif (opt==\'extBIC\') {\r\n+\t\te<--log10(ppoints(nrow(x$opt_extBIC$out)))\r\n+\t\to<--log10(sort(x$opt_extBIC$out$pval))\r\n+\t\tmaxp<-ceiling(max(o))\r\n+\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to extBIC\'))\r\n+\t\tabline(0,1,col="dark grey")}\r\n+\telse if (opt==\'mbonf\') {\r\n+\t\te<--log10(ppoints(nrow(x$opt_mbonf$out)))\r\n+\t\to<--log10(sort(x$opt_mbonf$out$pval))\r\n+\t\tmaxp<-ceiling(max(o))\r\n+\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to mbonf\'))\r\n+\t\tabline(0,1,col="dark grey")}\r\n+\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n+\r\n+\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 source_library/mlmm1.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/source_library/mlmm1.r Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| b"@@ -0,0 +1,486 @@\n+##############################################################################################################################################\r\n+###MLMM - Multi-Locus Mixed Model\r\n+###SET OF FUNCTIONS TO CARRY GWAS CORRECTING FOR POPULATION STRUCTURE WHILE INCLUDING COFACTORS THROUGH A STEPWISE-REGRESSION APPROACH\r\n+#######\r\n+#\r\n+##note: require EMMA\r\n+#library(emma)\r\n+#source('emma.r')\r\n+#\r\n+##REQUIRED DATA & FORMAT\r\n+#\r\n+#PHENOTYPE - Y: a vector of length m, with names(Y)=individual names\r\n+#GENOTYPE - X: a n by m matrix, where n=number of individuals, m=n umber of SNPs, with rownames(X)=individual names, and colnames(X)=SNP names\r\n+#KINSHIP - K: a n by n matrix, with rownames(K)=colnames(K)=individual names\r\n+#each of these data being sorted in the same way, according to the individual name\r\n+#\r\n+##FOR PLOTING THE GWAS RESULTS\r\n+#SNP INFORMATION - snp_info: a data frame having at least 3 columns:\r\n+# - 1 named 'SNP', with SNP names (same as colnames(X)),\r\n+# - 1 named 'Chr', with the chromosome number to which belong each SNP\r\n+# - 1 named 'Pos', with the position of the SNP onto the chromosome it belongs to.\r\n+#######\r\n+#\r\n+##FUNCTIONS USE\r\n+#save this file somewhere on your computer and source it!\r\n+#source('path/mlmm.r')\r\n+#\r\n+###FORWARD + BACKWARD ANALYSES\r\n+#mygwas<-mlmm(Y,X,K,nbchunks,maxsteps)\r\n+#X,Y,K as described above\r\n+#nbchunks: an integer defining the number of chunks of to run the analysis, allows to decrease the memory usage ==> minimum=2, increase it if you do not have enough memory\r\n+#maxsteps: maximum number of steps desired in the forward approach. The forward approach breaks automatically once the pseudo-heritability is close to 0,\r\n+#\t\t\thowever to avoid doing too many steps in case the pseudo-heritability does not reach a value close to 0, this parameter is also used.\r\n+#\t\t\tIt's value must be specified as an integer >= 3\r\n+#\r\n+###RESULTS\r\n+#\r\n+##STEPWISE TABLE\r\n+#mygwas$step_table\r\n+#\r\n+##PLOTS\r\n+#\r\n+##PLOTS FORM THE FORWARD TABLE\r\n+#plot_step_table(mygwas,type=c('h2','maxpval','BIC','extBIC'))\r\n+#\r\n+##RSS PLOT\r\n+#plot_step_RSS(mygwas)\r\n+#\r\n+##GWAS MANHATTAN PLOTS\r\n+#\r\n+#FORWARD STEPS\r\n+#plot_fwd_GWAS(mygwas,step,snp_info,pval_filt)\r\n+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#\r\n+#OPTIMAL MODELS\r\n+#Automatic identification of the optimal models within the forwrad-backward models according to the extendedBIC or multiple-bonferonni criteria\r\n+#\r\n+#plot_opt_GWAS(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#\r\n+##GWAS MANHATTAN PLOT ZOOMED IN A REGION OF INTEREST\r\n+#plot_fwd_region(mygwas,step,snp_info,pval_filt,chrom,pos1,pos2)\r\n+#step=the step to be plotted in the forward approach, where 1 is the EMMAX scan (no cofactor)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#chrom is an integer specifying the chromosome on which the region of interest is\r\n+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n+#\r\n+#plot_opt_region(mygwas,opt=c('extBIC','mbonf'),snp_info,pval_filt,chrom,pos1,pos2)\r\n+#snp_info as described above\r\n+#pval_filt=a p-value threshold for filtering the output, only p-vals below this threshold will be displayed in the plot\r\n+#chrom is an integer specifying the chromosome on which the region of interest is\r\n+#pos1, pos2 are integers delimiting the region of interest in the same unit as Pos in snp_info\r\n+#\r\n+##QQPLOTS of pvalues\r\n+#qqplot_fwd_GWAS(mygwas,nsteps)\r\n+#nsteps=maximum number of forward steps to be displayed\r\n+#\r\n+#qqplot_opt_GWAS(mygwas,opt=c('extBIC'"..b'output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7) \r\n+ legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))\r\n+ } else {\r\n+ legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))\r\n+ }\r\n+}\r\n+ \r\n+plot_region<-function(x,chrom,pos1,pos2){\r\n+\tregion<-subset(x$out,Chr==chrom & Pos>=pos1 & Pos <=pos2)\r\n+\tregion$col<- if (chrom %% 2 == 0) {\'gray60\'} else {\'gray10\'}\r\n+\tregion$col[which(region$SNP %in% x$cof)]<-\'red\'\r\n+\tplot(region$Pos,-log10(region$pval),type=\'p\',pch=20,main=paste(\'chromosome\',chrom,sep=\'\'),xlab=\'position (bp)\',ylab=\'-log10(pval)\',col=region$col,xlim=c(pos1,pos2))\r\n+\tabline(h=x$bonf_thresh,lty=3,col=\'black\')}\r\n+\r\n+\r\n+plot_fwd_GWAS<-function(x,step,snp_info,pval_filt) {\r\n+\tstopifnot(step<=length(x$pval_step))\r\n+\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n+\tplot_GWAS(output)}\r\n+\r\n+plot_fwd_region<-function(x,step,snp_info,pval_filt,chrom,pos1,pos2) {\r\n+\tstopifnot(step<=length(x$pval_step))\r\n+\toutput<-list(out=subset(merge(snp_info,x$pval_step[[step]]$out,by=\'SNP\'),pval<=pval_filt),cof=x$pval_step[[step]]$cof,bonf_thresh=x$bonf_thresh)\r\n+\tplot_region(output,chrom,pos1,pos2)}\r\n+\r\n+\r\n+plot_opt_GWAS<-function(x,opt,snp_info,pval_filt) {\r\n+\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_GWAS(output)}\r\n+\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_GWAS(output)}\r\n+\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n+\r\n+plot_opt_region<-function(x,opt,snp_info,pval_filt,chrom,pos1,pos2) {\r\n+\tif (opt==\'extBIC\') {output<-list(out=subset(merge(snp_info,x$opt_extBIC$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_extBIC$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_region(output,chrom,pos1,pos2)}\r\n+\telse if (opt==\'mbonf\') {output<-list(out=subset(merge(snp_info,x$opt_mbonf$out,by=\'SNP\'),pval<=pval_filt),cof=x$opt_mbonf$cof,bonf_thresh=x$bonf_thresh)\r\n+\t\tplot_region(output,chrom,pos1,pos2)}\r\n+\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n+\r\n+\r\n+qqplot_fwd_GWAS<-function(x,nsteps){\r\n+\tstopifnot(nsteps<=length(x$pval_step))\r\n+\te<--log10(ppoints(nrow(x$pval_step[[1]]$out)))\r\n+\tostep<-list()\r\n+\tostep[[1]]<--log10(sort(x$pval_step[[1]]$out$pval))\r\n+\tfor (i in 2:nsteps) {ostep[[i]]<--log10(sort(x$pval_step[[i]]$out$pval))}\r\n+\r\n+\tmaxp<-ceiling(max(unlist(ostep)))\r\n+\r\n+\tplot(e,ostep[[1]],type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp))\r\n+\tabline(0,1,col="dark grey")\r\n+\r\n+\tfor (i in 2:nsteps) {\r\n+\tpar(new=T)\r\n+\tplot(e,ostep[[i]],type=\'b\',pch=20,cex=0.8,col=i,axes=\'F\',xlab=\'\',ylab=\'\',xlim=c(0,max(e)+1),ylim=c(0,maxp))}\r\n+\tlegend(0,maxp,lty=1,pch=20,col=c(1:length(ostep)),paste(c(0:(length(ostep)-1)),\'cof\',sep=\' \'))\r\n+}\r\n+\r\n+qqplot_opt_GWAS<-function(x,opt){\r\n+\tif (opt==\'extBIC\') {\r\n+\t\te<--log10(ppoints(nrow(x$opt_extBIC$out)))\r\n+\t\to<--log10(sort(x$opt_extBIC$out$pval))\r\n+\t\tmaxp<-ceiling(max(o))\r\n+\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to extBIC\'))\r\n+\t\tabline(0,1,col="dark grey")}\r\n+\telse if (opt==\'mbonf\') {\r\n+\t\te<--log10(ppoints(nrow(x$opt_mbonf$out)))\r\n+\t\to<--log10(sort(x$opt_mbonf$out$pval))\r\n+\t\tmaxp<-ceiling(max(o))\r\n+\t\tplot(e,o,type=\'b\',pch=20,cex=0.8,col=1,xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,max(e)+1),ylim=c(0,maxp),main=paste(\'optimal model according to mbonf\'))\r\n+\t\tabline(0,1,col="dark grey")}\r\n+\telse {cat(\'error! \\n opt must be extBIC or mbonf\')}}\r\n+\r\n+\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 source_library/plot_MLMM.Ago.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/source_library/plot_MLMM.Ago.r Mon Apr 16 08:50:05 2018 -0400 |
| [ |
| @@ -0,0 +1,36 @@ +## fonction ########################################################################### +plot_MLMM.Ago<-function(x) { + output1 <-x$out[order(x$out$Pos),] + output_ok<-output1[order(output1$Chr),] + + maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200))) + plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2)) +# plot_col<-c('blue','darkgreen','red','cyan','purple') + size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x + a<-rep(maxpos[1],size[1]) + b<-rep(plot_col[1],size[1]) + for (i in 2:max(unique(output_ok$Chr))){ + a<-c(a,rep(maxpos[i],size[i])) + b<-c(b,rep(plot_col[i],size[i]))} + + output_ok$xpos<-output_ok$Pos+a + output_ok$col<-b + output_ok$col[output_ok$mk=='qtl']<-'cyan' + output_ok$col[output_ok$SNP %in% x$cof]<-'red' + + d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2 + + plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome') + axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr)))) + xline(output_ok$xpos[output_ok$mk=='qtl'], col='cyan', lwd=0.1) + abline(h=x$bonf_thresh,lty=3,col='black') + + + if (length(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]) > 0) { + text(output_ok$xpos[-log10(output_ok$pval) > x$bonf_thresh], -log10(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7) + legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" ")) + } else { + legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" ")) + } + +} |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/PCs.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/PCs.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| b'@@ -0,0 +1,501 @@\n+Ind_id\tPC1\tPC2\tPC3\tPC4\tPC5\tPC6\tPC7\tPC8\tPC9\tPC10\r\n+Ind1\t-178.865919402083\t52.7627738156725\t-141.748945536388\t32.058224360871\t17.5543496781891\t-14.3548988388948\t-10.2007596874656\t-27.5522180438276\t53.8361897205006\t-18.6475470130354\r\n+Ind2\t17.2004199310976\t-35.0481518976658\t2.99233867314198\t4.81139913582997\t10.2228963692622\t-20.3297331439499\t39.2713692748573\t24.7110877292016\t-24.875136340639\t-79.7254551692405\r\n+Ind3\t-0.134476678588869\t-40.4809873117503\t-21.6926836525902\t-60.4011699911492\t-22.0440390400408\t79.7396191270007\t-5.03699286956945\t1.44260806327655\t-18.5673225023224\t0.783908106338091\r\n+Ind4\t59.6311928105328\t2.40615133725044\t-8.46928822591139\t4.73119626189056\t-57.1868128931894\t-20.7915525617803\t30.7364511052182\t47.9466786336674\t2.06990342944568\t-18.0535676719438\r\n+Ind5\t103.891300912847\t21.8037845187655\t-30.6709667453021\t1.37837230460589\t15.0719778307841\t8.38159537670945\t83.5520948354437\t70.186244359573\t-74.0537682988908\t-186.781952769146\r\n+Ind6\t109.589582004282\t-118.667027609925\t-56.5722087552307\t48.7459288483199\t43.1571558829337\t7.79129644694251\t-17.772465642767\t-18.1684020333467\t6.43049062690154\t18.2609441658189\r\n+Ind7\t105.75300681181\t-93.7404190707278\t-51.0560736617596\t45.3952340632916\t8.61940092138935\t-3.32544142648814\t8.76332422786606\t17.82243564091\t5.74090720937479\t-5.4054735161598\r\n+Ind8\t159.117919241928\t210.416426572268\t4.20702925818152\t-11.8376448321144\t53.2323518155232\t10.3499980398735\t-15.1480685184512\t1.00301961757128\t17.2748848820369\t23.6530269862681\r\n+Ind9\t-50.1469406712458\t-38.2756828481621\t71.5460249979154\t-76.8836948044296\t41.1850635495537\t-49.8829649979189\t53.0762145097006\t-44.8704271004554\t-34.5416781774017\t21.333991327445\r\n+Ind10\t-23.6824581856269\t23.4679655693919\t25.5094640956877\t-9.2821017205585\t-36.0153179954092\t-6.09111990033547\t37.4935851537382\t7.31118681274309\t12.8801360250816\t-21.1778668967879\r\n+Ind11\t-70.5695919961539\t-17.9167749647058\t63.0214219947379\t5.39107768462564\t20.2888030311226\t-9.45803953937197\t8.87556183720302\t-19.4761604377054\t32.2587991979664\t-25.0662549787725\r\n+Ind12\t77.2077765285411\t-62.659578207814\t-27.8657805489511\t20.7987141108637\t-1.41821017110975\t7.5993766669185\t17.9690999522708\t22.9642475227629\t-15.3421596468532\t-39.1689172788894\r\n+Ind13\t-67.3091822321392\t2.66980623558726\t80.044052393184\t24.6724122096384\t-7.48210168827798\t-4.46029731717329\t7.75547732175856\t-14.0675127206306\t32.8839994806595\t-10.4795635690372\r\n+Ind14\t-184.034884860198\t56.7224761666409\t-67.474111495687\t66.9996283765167\t5.47192566420099\t3.89211715728053\t-27.6339252234967\t-31.1198732901585\t136.24294045799\t-44.7542900589228\r\n+Ind15\t105.428042508281\t-113.589279126332\t-51.8683518173688\t57.9619804375842\t37.8043886385277\t-9.45565351031311\t-19.8829717462431\t-30.0885630584139\t16.4261500015328\t18.7011808636841\r\n+Ind16\t-166.999615746489\t43.6671290195502\t-60.2943993572178\t46.0254896111197\t25.3231293734176\t-16.6610488229775\t-20.8346628454123\t-6.08489275545223\t14.8219858951933\t-31.9509911024227\r\n+Ind17\t-56.3817531071442\t-59.4028530568672\t78.9406978235452\t-107.360391080384\t67.8224855343535\t-46.7656143929592\t-106.555590573541\t87.1230040934812\t24.4028645049117\t-13.1509389429015\r\n+Ind18\t-39.7809078624736\t-25.8176912681573\t73.4765098821242\t-22.7023360738786\t14.6723870006094\t-21.6956525811808\t23.8367824362572\t-15.2260511829874\t23.4711883636939\t-14.7132705886433\r\n+Ind19\t-78.5680239119424\t18.2235502556336\t74.7093453039448\t34.4925314908252\t-6.56036047030315\t-10.2793681564071\t22.4651313808755\t-13.881271526365\t36.0293206033664\t-14.7533674484279\r\n+Ind20\t48.2470569421969\t-3.69020553055465\t-9.59636579892298\t-6.69680116703723\t-29.8614802157303\t-9.34736756131693\t32.3419807777327\t32.493014414624\t2.23994486173309\t-13.5480000325967\r\n+Ind21\t-60.6918386867323\t-46.2391580805452\t82.3767657175553\t-97.8344788950968\t42.5792486239748\t-46.1392500754955\t36.3208810132587\t-74.4779090496684\t-36.8763303734952\t19.3809530214196\r\n+Ind22\t-20.1597864683843\t-7.23057579635946\t53.6695295014678\t-26.3081178695212\t14.3467873437096\t-34.1297211748669'..b'\t52.5665305760732\t-65.3330748667134\t42.5311878441339\t-24.464174943151\t-3.07473433090338\t0.554397998009222\t-5.28165770901293\t-2.78870111837357\r\n+Ind480\t-4.04752303092518\t-73.6992501575367\t-65.4031344401346\t-164.801942671165\t-104.089690224183\t396.297584404712\t-66.4385559351399\t-99.0230888742127\t38.4508361256368\t-3.97361256493684\r\n+Ind481\t61.5862811316281\t-78.0786733152836\t-15.3296069848823\t-13.6289850803238\t5.14773569036415\t24.2749808058283\t-8.78156065647312\t-2.87527282007219\t-12.9592131231043\t-3.11169624125549\r\n+Ind482\t-60.408002836006\t-37.7904990541033\t75.500818076423\t-93.2482190610767\t42.3896796098016\t-31.2928678476235\t38.3361659058781\t-65.7777184353998\t-31.8359086971969\t26.7379688844134\r\n+Ind483\t111.199202280542\t-124.80420233033\t-60.8380227039332\t64.8441823398563\t54.8715109960327\t-25.4010373310897\t-12.9241890463721\t-18.980657154477\t15.8286980189606\t29.0780788349656\r\n+Ind484\t-180.948979506013\t62.6368939791308\t-214.916990365476\t-2.75712119030329\t16.3524253928772\t-32.8245275261777\t-5.24395924099511\t0.917462670165947\t-14.4836038771732\t28.009930300943\r\n+Ind485\t-48.9586054213367\t-40.4437207405978\t79.5207365630228\t-84.1510490182313\t42.3841568298077\t-57.8341012798792\t59.7480876119257\t-66.9341128043664\t-42.3657697351523\t14.775299966591\r\n+Ind486\t-159.666890251042\t31.4966506310633\t-139.131632905916\t-11.2618218296909\t14.7453387872702\t-10.205452714936\t-12.124864689413\t18.1899241017654\t-8.72925727860455\t24.019668975405\r\n+Ind487\t50.5880802752021\t48.0170747575094\t6.62029768869207\t6.49587626282199\t-108.205391743711\t-50.4660828797446\t-18.4304604339707\t-1.92138539408283\t-4.40803505129628\t7.56963246091155\r\n+Ind488\t-208.657797405578\t83.4414645106006\t-270.60556716331\t-9.22849863470282\t19.4750154595157\t-40.0180837031553\t1.69169862979747\t2.59575037957937\t-30.9897831274821\t28.5650087815123\r\n+Ind489\t109.008174636546\t-127.145014483296\t-53.7208205104242\t59.000400588283\t52.5582108214735\t-13.4425977641921\t-23.0408546584087\t-32.3408372189686\t19.4529599178797\t32.4093332909555\r\n+Ind490\t120.716268506823\t-121.623287466096\t-60.2563278213462\t71.1137507086527\t52.5730997111148\t-15.4377355868701\t-25.4373797449979\t-30.0019237706319\t19.1811467119508\t36.5431608680519\r\n+Ind491\t53.4116590621561\t-47.4233860769833\t-17.2170658213648\t3.53813216420503\t15.2304866328875\t8.67380978028018\t26.3321125878518\t10.5013916619316\t-21.3298007515567\t-23.7400513114709\r\n+Ind492\t-53.1664781930866\t-70.581036023938\t93.645992497218\t-151.590091323412\t95.8363803040874\t-49.2920792895683\t-192.299333969135\t155.536237355319\t31.3435787246676\t-20.6898528099124\r\n+Ind493\t-5.19680901320142\t-14.0934775822778\t-7.95320505870183\t-37.6033484002868\t-65.3243315887459\t33.5845682283006\t40.2595466161277\t25.0513283224683\t-4.173645727916\t14.9527489104772\r\n+Ind494\t-97.0678356035457\t14.3683418908876\t106.972882512285\t117.433956227568\t-2.05375973333015\t52.0427486776749\t-29.1707532151143\t20.7121723820673\t-69.7977680085945\t15.0019034321513\r\n+Ind495\t20.2280926278577\t26.3525021428749\t-7.61248590214303\t-34.8739039972602\t-61.9484758979425\t44.9434884827441\t47.3618447482387\t28.8415799262384\t0.572917654609841\t4.6821282141236\r\n+Ind496\t-87.5401592241228\t-3.8746585095172\t-114.784466762102\t-34.8106027846935\t39.9049599166775\t-17.5510387083759\t12.0478069658709\t-0.826997872575238\t-69.122613128298\t2.31363417720061\r\n+Ind497\t-9.58901759291861\t-52.3711752519947\t-34.0558865517471\t-54.8315082991333\t-28.2563953822263\t74.2474833700731\t-1.6250241765911\t17.3448111549602\t2.71321560989128\t4.8116619479419\r\n+Ind498\t102.522038959597\t-117.460334015358\t-52.9498836241136\t66.5117949960196\t40.7147989150976\t-14.4738432511714\t-7.03603829584866\t-9.75651513912602\t9.48252046472213\t7.01133398260746\r\n+Ind499\t174.49863569711\t267.625756275349\t7.01439154444572\t-17.081552800905\t100.001054773953\t31.9244766956512\t-50.0344497868791\t-53.0688352612071\t-20.1600763500728\t-21.4894632239011\r\n+Ind500\t-37.6409993591417\t-40.6791977977316\t57.7896893905105\t-60.1120460885041\t15.9332832089983\t-1.48808724253164\t11.3154567784946\t-24.2489700785127\t-19.2823043816803\t7.88319137420756\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/genot.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/genot.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| b'@@ -0,0 +1,501 @@\n+Ind_id\tSNP1\tSNP2\tSNP3\tSNP4\tSNP5\tSNP6\tSNP7\tSNP8\tSNP9\tSNP10\tSNP11\tSNP12\tSNP13\tSNP14\tSNP15\tSNP16\tSNP17\tSNP18\tSNP19\tSNP20\tSNP21\tSNP22\tSNP23\tSNP24\tSNP25\tSNP26\tSNP27\tSNP28\tSNP29\tSNP30\tSNP31\tSNP32\tSNP33\tSNP34\tSNP35\tSNP36\tSNP37\tSNP38\tSNP39\tSNP40\tSNP41\tSNP42\tSNP43\tSNP44\tSNP45\tSNP46\tSNP47\tSNP48\tSNP49\tSNP50\tSNP51\tSNP52\tSNP53\tSNP54\tSNP55\tSNP56\tSNP57\tSNP58\tSNP59\tSNP60\tSNP61\tSNP62\tSNP63\tSNP64\tSNP65\tSNP66\tSNP67\tSNP68\tSNP69\tSNP70\tSNP71\tSNP72\tSNP73\tSNP74\tSNP75\tSNP76\tSNP77\tSNP78\tSNP79\tSNP80\tSNP81\tSNP82\tSNP83\tSNP84\tSNP85\tSNP86\tSNP87\tSNP88\tSNP89\tSNP90\tSNP91\tSNP92\tSNP93\tSNP94\tSNP95\tSNP96\tSNP97\tSNP98\tSNP99\tSNP100\tSNP101\tSNP102\tSNP103\tSNP104\tSNP105\tSNP106\tSNP107\tSNP108\tSNP109\tSNP110\tSNP111\tSNP112\tSNP113\tSNP114\tSNP115\tSNP116\tSNP117\tSNP118\tSNP119\tSNP120\tSNP121\tSNP122\tSNP123\tSNP124\tSNP125\tSNP126\tSNP127\tSNP128\tSNP129\tSNP130\tSNP131\tSNP132\tSNP133\tSNP134\tSNP135\tSNP136\tSNP137\tSNP138\tSNP139\tSNP140\tSNP141\tSNP142\tSNP143\tSNP144\tSNP145\tSNP146\tSNP147\tSNP148\tSNP149\tSNP150\tSNP151\tSNP152\tSNP153\tSNP154\tSNP155\tSNP156\tSNP157\tSNP158\tSNP159\tSNP160\tSNP161\tSNP162\tSNP163\tSNP164\tSNP165\tSNP166\tSNP167\tSNP168\tSNP169\tSNP170\tSNP171\tSNP172\tSNP173\tSNP174\tSNP175\tSNP176\tSNP177\tSNP178\tSNP179\tSNP180\tSNP181\tSNP182\tSNP183\tSNP184\tSNP185\tSNP186\tSNP187\tSNP188\tSNP189\tSNP190\tSNP191\tSNP192\tSNP193\tSNP194\tSNP195\tSNP196\tSNP197\tSNP198\tSNP199\tSNP200\tSNP201\tSNP202\tSNP203\tSNP204\tSNP205\tSNP206\tSNP207\tSNP208\tSNP209\tSNP210\tSNP211\tSNP212\tSNP213\tSNP214\tSNP215\tSNP216\tSNP217\tSNP218\tSNP219\tSNP220\tSNP221\tSNP222\tSNP223\tSNP224\tSNP225\tSNP226\tSNP227\tSNP228\tSNP229\tSNP230\tSNP231\tSNP232\tSNP233\tSNP234\tSNP235\tSNP236\tSNP237\tSNP238\tSNP239\tSNP240\tSNP241\tSNP242\tSNP243\tSNP244\tSNP245\tSNP246\tSNP247\tSNP248\tSNP249\tSNP250\tSNP251\tSNP252\tSNP253\tSNP254\tSNP255\tSNP256\tSNP257\tSNP258\tSNP259\tSNP260\tSNP261\tSNP262\tSNP263\tSNP264\tSNP265\tSNP266\tSNP267\tSNP268\tSNP269\tSNP270\tSNP271\tSNP272\tSNP273\tSNP274\tSNP275\tSNP276\tSNP277\tSNP278\tSNP279\tSNP280\tSNP281\tSNP282\tSNP283\tSNP284\tSNP285\tSNP286\tSNP287\tSNP288\tSNP289\tSNP290\tSNP291\tSNP292\tSNP293\tSNP294\tSNP295\tSNP296\tSNP297\tSNP298\tSNP299\tSNP300\tSNP301\tSNP302\tSNP303\tSNP304\tSNP305\tSNP306\tSNP307\tSNP308\tSNP309\tSNP310\tSNP311\tSNP312\tSNP313\tSNP314\tSNP315\tSNP316\tSNP317\tSNP318\tSNP319\tSNP320\tSNP321\tSNP322\tSNP323\tSNP324\tSNP325\tSNP326\tSNP327\tSNP328\tSNP329\tSNP330\tSNP331\tSNP332\tSNP333\tSNP334\tSNP335\tSNP336\tSNP337\tSNP338\tSNP339\tSNP340\tSNP341\tSNP342\tSNP343\tSNP344\tSNP345\tSNP346\tSNP347\tSNP348\tSNP349\tSNP350\tSNP351\tSNP352\tSNP353\tSNP354\tSNP355\tSNP356\tSNP357\tSNP358\tSNP359\tSNP360\tSNP361\tSNP362\tSNP363\tSNP364\tSNP365\tSNP366\tSNP367\tSNP368\tSNP369\tSNP370\tSNP371\tSNP372\tSNP373\tSNP374\tSNP375\tSNP376\tSNP377\tSNP378\tSNP379\tSNP380\tSNP381\tSNP382\tSNP383\tSNP384\tSNP385\tSNP386\tSNP387\tSNP388\tSNP389\tSNP390\tSNP391\tSNP392\tSNP393\tSNP394\tSNP395\tSNP396\tSNP397\tSNP398\tSNP399\tSNP400\tSNP401\tSNP402\tSNP403\tSNP404\tSNP405\tSNP406\tSNP407\tSNP408\tSNP409\tSNP410\tSNP411\tSNP412\tSNP413\tSNP414\tSNP415\tSNP416\tSNP417\tSNP418\tSNP419\tSNP420\tSNP421\tSNP422\tSNP423\tSNP424\tSNP425\tSNP426\tSNP427\tSNP428\tSNP429\tSNP430\tSNP431\tSNP432\tSNP433\tSNP434\tSNP435\tSNP436\tSNP437\tSNP438\tSNP439\tSNP440\tSNP441\tSNP442\tSNP443\tSNP444\tSNP445\tSNP446\tSNP447\tSNP448\tSNP449\tSNP450\tSNP451\tSNP452\tSNP453\tSNP454\tSNP455\tSNP456\tSNP457\tSNP458\tSNP459\tSNP460\tSNP461\tSNP462\tSNP463\tSNP464\tSNP465\tSNP466\tSNP467\tSNP468\tSNP469\tSNP470\tSNP471\tSNP472\tSNP473\tSNP474\tSNP475\tSNP476\tSNP477\tSNP478\tSNP479\tSNP480\tSNP481\tSNP482\tSNP483\tSNP484\tSNP485\tSNP486\tSNP487\tSNP488\tSNP489\tSNP490\tSNP491\tSNP492\tSNP493\tSNP494\tSNP495\tSNP496\tSNP497\tSNP498\tSNP499\tSNP500\tSNP501\tSNP502\tSNP503\tSNP504\tSNP505\tSNP506\tSNP507\tSNP508\tSNP509\tSNP510\tSNP511\tSNP512\tSNP513\tSNP514\tSNP515\tSNP516\tSNP517\tSNP518\tSNP519\tSNP520\tSNP521\tSNP522\tSNP523\tSNP524\tSNP525\tSNP526\tSNP527\tSNP528\tSNP529\tSNP530\tSNP531\tSNP532\tSNP533\tSNP534\tSNP535\tSNP536\tSNP537\tSNP538\tSNP539\tSNP540\tSNP541\tSNP542\tSNP543\tSNP544\tSNP545\tSNP546\tSNP547\tSNP548\tSNP549\tSNP550\tSNP551\tSNP552\tSNP553\tSNP554\tSNP555\tSNP556\tSNP557\tSNP558\tSNP559\tSNP560\tSNP561\tSNP562\tSNP563\tSNP564\tSNP565\tSNP566\tSNP567\tSNP568\tSNP569\tSNP570\tSNP571\tSNP572\tSNP573\tSNP574\tSNP575\tSNP576\tSNP577\tSNP578\tSNP579\tSNP580\tSNP581\tSNP582\tSNP583\t'..b'\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t2\t0\t2\t0\tNA\t0\t0\t2\t0\t0\t2\t0\t2\t0\t2\t0\t2\t2\t0\t0\t2\t2\t2\t0\t2\t0\t2\t0\t2\t2\t0\t0\t0\t2\t2\t0\t2\t0\t0\t0\t0\t2\t2\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t0\t2\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t2\t2\t0\t2\t0\t2\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t2\t2\t0\t0\t2\t2\tNA\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t2\t2\t0\t2\t0\t2\t0\t0\t2\t0\t2\t0\t2\t2\t2\t2\t2\t2\t0\t0\t0\t2\t2\t2\t0\t0\t0\t2\t0\t0\t0\t2\t2\t0\t0\t2\t2\t0\t0\t2\t0\t0\t0\t2\t2\t0\t0\t2\t2\t0\t2\t2\t0\t0\t0\t2\t2\tNA\t0\t0\t0\t0\t0\t0\t0\t2\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t2\t0\t2\t2\t2\t0\t2\t2\t0\t2\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t2\t0\t0\t2\t2\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t0\t2\t2\t2\t0\t0\t0\t2\t2\t0\t0\t2\t0\t0\t0\t0\t2\t2\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t2\t2\t0\t0\t0\t2\t0\t0\t0\t0\t2\t0\t2\t0\t2\t2\t2\t2\t0\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t2\t0\t2\t0\t2\t2\t2\t0\t2\t2\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t2\t0\t0\t2\t0\t0\t0\t2\t2\t0\t0\t0\t0\t2\t2\t0\t2\t2\t2\t2\t2\t0\t0\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t2\t2\t2\t2\t0\t2\t2\t0\t2\t0\t2\t0\t0\t0\t2\t0\t2\t0\t2\t0\tNA\t0\t0\t0\t2\t0\t0\t0\t2\t2\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t0\t2\t2\t0\t2\t0\t0\t2\t2\t0\t0\t0\t0\t0\t2\t2\t2\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t2\t2\t0\t0\t0\t2\t2\t2\t2\t2\t0\t2\t0\t2\t2\t0\t0\t0\t2\t0\t2\t2\t2\t0\t0\tNA\t0\t0\t0\t0\t0\t0\t2\t2\t0\t2\t2\t2\t2\t2\t0\t2\t2\t0\t2\t0\t0\t2\t2\t0\t0\t2\t0\t0\t2\t0\t0\t0\t2\t0\t2\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t0\t2\t2\t0\t2\t0\t2\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t2\t2\t0\t2\t0\t2\t2\t0\t0\tNA\t0\t0\t0\t0\t0\t2\t2\t2\t2\t2\t2\t0\t0\t2\t0\t0\t0\t2\t0\t2\t0\t0\t0\t2\t2\t2\t0\t2\t2\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t2\t2\t2\t2\t0\t0\t2\t2\t0\t2\t0\t2\t2\t0\t0\t0\t0\t0\t2\t0\t2\t2\t2\t2\t0\t2\t2\t2\t2\t2\t0\t2\t0\t2\t2\t2\t0\t0\t2\t0\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t0\t0\t0\t2\t2\t0\t0\t2\t2\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t2\t0\t2\t2\t0\t2\t0\t0\t2\t0\t0\t0\t2\t2\t0\t2\t0\t0\tNA\t0\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t2\t0\t2\t2\t0\t2\t0\t0\t2\t0\t0\t2\t2\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t2\t2\t2\t2\t2\t0\t2\t0\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t2\t2\t2\t2\t2\t0\t2\t0\t2\t0\t2\t2\t0\t2\t2\t0\t2\t2\t0\t0\t0\t2\t0\t2\t0\t0\t0\t2\t2\t2\t0\t2\t0\t0\t0\t2\t0\t2\t0\t2\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t0\t2\t0\t2\t0\t0\t2\tNA\t2\t0\t2\t2\t2\t0\t2\t2\t0\t0\t0\t2\t0\t0\t2\t0\t0\t2\t0\t0\t0\t2\t2\t0\t2\tNA\t0\t0\t2\t2\t0\t0\t0\t0\t2\t0\t2\t2\t0\t0\t0\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t2\t2\t2\t2\t0\t0\t2\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t2\t0\t2\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\tNA\t2\t0\t2\t2\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t0\t0\t2\t2\t0\t2\t0\t2\t2\t0\t2\t0\t0\t2\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t2\t2\t0\t0\t0\t2\t0\t0\t2\tNA\t0\t0\t0\t2\t2\t0\t0\t0\t0\t0\t2\t2\t0\t0\t0\t0\t2\t0\t0\t0\t2\t2\t2\t0\t0\t0\t0\t2\t2\t0\t2\t2\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t2\t2\t2\t2\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t2\t0\t2\t2\t0\t0\t2\t2\t0\t2\t2\t2\t2\t0\t2\t0\t0\t0\t2\t2\t0\t0\t0\t2\t0\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t2\t2\t2\t0\t0\t2\t2\t2\t0\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t2\t0\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t2\t2\t0\t0\t2\t0\t0\t0\t0\t0\t2\t2\t2\t0\t0\t2\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t2\t0\t0\t0\t2\t0\t0\t2\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t2\t0\t0\t0\t2\t2\t2\t0\t2\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t2\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t2\t2\t2\t2\t2\t0\t0\t2\t2\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t2\t2\t0\t0\t2\t2\t0\t0\t0\t2\t0\t2\t0\t0\t2\t2\t0\t0\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t2\t2\t0\t0\t2\t2\t0\t2\t0\t0\t2\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t2\t0\t2\t0\t2\t2\t0\t2\t0\t2\t0\t2\t0\t0\t2\t0\t0\t2\t0\t2\t0\t0\t0\tNA\t2\t2\t2\t2\t0\t2\t0\t0\t2\t2\t0\t2\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t2\t0\t2\t2\t0\t0\t0\t2\t2\t2\t0\t0\t0\t0\t2\t0\t2\t2\t2\t2\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t2\t2\t2\t2\t0\t2\t0\t0\t0\t2\t2\t0\t0\t2\t0\t2\t0\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t2\t0\t0\t0\t2\t2\t0\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t0\t2\t0\t0\t2\t2\t0\t0\t2\t2\t2\t2\t2\t0\t2\t0\t0\t2\t0\t0\t2\t2\t0\t2\t2\t0\t2\t0\t0\t0\t0\t2\t0\t2\t2\t2\t2\t2\t0\t0\t2\t2\t2\t2\t0\t2\t0\t0\t2\t2\t0\t2\t2\t0\t0\t0\t2\t0\t2\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t2\t0\t0\t2\t0\t0\t0\t2\t2\t2\t0\t0\t0\t0\t2\t0\tNA\t0\t2\t0\t2\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t0\t2\t2\t2\t0\t0\t0\t0\t2\t0\t0\t0\t2\t2\t2\t0\t0\t2\t0\t2\t0\t0\t2\t2\t2\t2\t2\t0\t0\t2\t2\t0\t0\t0\t0\t2\t0\t0\t0\t2\t0\t2\t0\t0\t2\t2\t2\t0\t2\t2\t0\t2\t0\t2\t0\t2\t0\t0\t2\t2\t0\t2\t0\t2\t0\t2\t0\t2\t0\t0\t0\t2\t2\t0\t2\t2\t0\t0\t0\t2\t2\t0\t0\t0\t2\t0\t0\t2\t0\tNA\t2\t0\t2\t2\t2\t0\t2\t2\t2\t2\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t2\t0\t0\t2\t0\t0\t0\t0\t2\t0\t2\t2\t2\t0\t2\t0\t2\t0\t2\t0\t0\t2\t2\t2\tNA\t2\t2\t0\t2\t0\t2\t2\t0\t0\t2\t0\t0\t2\t2\t2\t0\t0\t0\t2\t0\tNA\t0\t0\t0\t2\t0\t0\t0\t0\t0\t2\t2\t0\t2\t0\t2\t2\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t2\t0\t2\t0\t0\t0\t2\t0\t0\t0\t0\t0\t0\t2\t0\t2\t0\t0\t2\t0\t0\t2\t0\t2\t2\t2\t0\t0\t0\t0\t0\t2\t2\t2\t0\t0\t2\t0\t0\t0\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t0\t2\t0\t0\t2\t2\t0\t0\t0\t2\t0\t2\t2\t0\t2\t2\t0\t0\t0\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t0\t0\t2\t2\t2\t0\t0\t2\t0\t0\t0\t0\t2\t2\t2\t0\t0\t2\t0\t0\t2\t2\t2\t0\t0\t2\t0\t0\t2\t0\t0\t2\t2\t2\t0\t2\t0\t0\t0\t2\t0\t0\t2\t0\t2\t0\t0\t0\t0\t0\t0\t2\t2\t0\t2\t0\t0\t0\t2\t2\t0\t2\t0\t0\t0\t0\t2\t2\t2\t2\t2\t0\t0\t0\t0\t0\t0\t0\t2\t2\t0\t0\t2\t0\t0\t2\t0\t2\t0\t2\t0\t0\t0\t0\n\\ No newline at end of file\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/kinship.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/kinship.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| b'@@ -0,0 +1,501 @@\n+Ind1\tInd2\tInd3\tInd4\tInd5\tInd6\tInd7\tInd8\tInd9\tInd10\tInd11\tInd12\tInd13\tInd14\tInd15\tInd16\tInd17\tInd18\tInd19\tInd20\tInd21\tInd22\tInd23\tInd24\tInd25\tInd26\tInd27\tInd28\tInd29\tInd30\tInd31\tInd32\tInd33\tInd34\tInd35\tInd36\tInd37\tInd38\tInd39\tInd40\tInd41\tInd42\tInd43\tInd44\tInd45\tInd46\tInd47\tInd48\tInd49\tInd50\tInd51\tInd52\tInd53\tInd54\tInd55\tInd56\tInd57\tInd58\tInd59\tInd60\tInd61\tInd62\tInd63\tInd64\tInd65\tInd66\tInd67\tInd68\tInd69\tInd70\tInd71\tInd72\tInd73\tInd74\tInd75\tInd76\tInd77\tInd78\tInd79\tInd80\tInd81\tInd82\tInd83\tInd84\tInd85\tInd86\tInd87\tInd88\tInd89\tInd90\tInd91\tInd92\tInd93\tInd94\tInd95\tInd96\tInd97\tInd98\tInd99\tInd100\tInd101\tInd102\tInd103\tInd104\tInd105\tInd106\tInd107\tInd108\tInd109\tInd110\tInd111\tInd112\tInd113\tInd114\tInd115\tInd116\tInd117\tInd118\tInd119\tInd120\tInd121\tInd122\tInd123\tInd124\tInd125\tInd126\tInd127\tInd128\tInd129\tInd130\tInd131\tInd132\tInd133\tInd134\tInd135\tInd136\tInd137\tInd138\tInd139\tInd140\tInd141\tInd142\tInd143\tInd144\tInd145\tInd146\tInd147\tInd148\tInd149\tInd150\tInd151\tInd152\tInd153\tInd154\tInd155\tInd156\tInd157\tInd158\tInd159\tInd160\tInd161\tInd162\tInd163\tInd164\tInd165\tInd166\tInd167\tInd168\tInd169\tInd170\tInd171\tInd172\tInd173\tInd174\tInd175\tInd176\tInd177\tInd178\tInd179\tInd180\tInd181\tInd182\tInd183\tInd184\tInd185\tInd186\tInd187\tInd188\tInd189\tInd190\tInd191\tInd192\tInd193\tInd194\tInd195\tInd196\tInd197\tInd198\tInd199\tInd200\tInd201\tInd202\tInd203\tInd204\tInd205\tInd206\tInd207\tInd208\tInd209\tInd210\tInd211\tInd212\tInd213\tInd214\tInd215\tInd216\tInd217\tInd218\tInd219\tInd220\tInd221\tInd222\tInd223\tInd224\tInd225\tInd226\tInd227\tInd228\tInd229\tInd230\tInd231\tInd232\tInd233\tInd234\tInd235\tInd236\tInd237\tInd238\tInd239\tInd240\tInd241\tInd242\tInd243\tInd244\tInd245\tInd246\tInd247\tInd248\tInd249\tInd250\tInd251\tInd252\tInd253\tInd254\tInd255\tInd256\tInd257\tInd258\tInd259\tInd260\tInd261\tInd262\tInd263\tInd264\tInd265\tInd266\tInd267\tInd268\tInd269\tInd270\tInd271\tInd272\tInd273\tInd274\tInd275\tInd276\tInd277\tInd278\tInd279\tInd280\tInd281\tInd282\tInd283\tInd284\tInd285\tInd286\tInd287\tInd288\tInd289\tInd290\tInd291\tInd292\tInd293\tInd294\tInd295\tInd296\tInd297\tInd298\tInd299\tInd300\tInd301\tInd302\tInd303\tInd304\tInd305\tInd306\tInd307\tInd308\tInd309\tInd310\tInd311\tInd312\tInd313\tInd314\tInd315\tInd316\tInd317\tInd318\tInd319\tInd320\tInd321\tInd322\tInd323\tInd324\tInd325\tInd326\tInd327\tInd328\tInd329\tInd330\tInd331\tInd332\tInd333\tInd334\tInd335\tInd336\tInd337\tInd338\tInd339\tInd340\tInd341\tInd342\tInd343\tInd344\tInd345\tInd346\tInd347\tInd348\tInd349\tInd350\tInd351\tInd352\tInd353\tInd354\tInd355\tInd356\tInd357\tInd358\tInd359\tInd360\tInd361\tInd362\tInd363\tInd364\tInd365\tInd366\tInd367\tInd368\tInd369\tInd370\tInd371\tInd372\tInd373\tInd374\tInd375\tInd376\tInd377\tInd378\tInd379\tInd380\tInd381\tInd382\tInd383\tInd384\tInd385\tInd386\tInd387\tInd388\tInd389\tInd390\tInd391\tInd392\tInd393\tInd394\tInd395\tInd396\tInd397\tInd398\tInd399\tInd400\tInd401\tInd402\tInd403\tInd404\tInd405\tInd406\tInd407\tInd408\tInd409\tInd410\tInd411\tInd412\tInd413\tInd414\tInd415\tInd416\tInd417\tInd418\tInd419\tInd420\tInd421\tInd422\tInd423\tInd424\tInd425\tInd426\tInd427\tInd428\tInd429\tInd430\tInd431\tInd432\tInd433\tInd434\tInd435\tInd436\tInd437\tInd438\tInd439\tInd440\tInd441\tInd442\tInd443\tInd444\tInd445\tInd446\tInd447\tInd448\tInd449\tInd450\tInd451\tInd452\tInd453\tInd454\tInd455\tInd456\tInd457\tInd458\tInd459\tInd460\tInd461\tInd462\tInd463\tInd464\tInd465\tInd466\tInd467\tInd468\tInd469\tInd470\tInd471\tInd472\tInd473\tInd474\tInd475\tInd476\tInd477\tInd478\tInd479\tInd480\tInd481\tInd482\tInd483\tInd484\tInd485\tInd486\tInd487\tInd488\tInd489\tInd490\tInd491\tInd492\tInd493\tInd494\tInd495\tInd496\tInd497\tInd498\tInd499\tInd500\n+Ind1\t0.980555602735708\t-0.0329727937328273\t-0.0262487415186466\t-0.0604104840755713\t-0.0684357122797762\t-0.078471107107226\t-0.0706737072888241\t-0.0755481118857582\t-0.012703922999043\t0.00235095437161519\t0.0673760868098489\t-0.0727124982715698\t-0.00272547502665574\t0.367287948459605\t-0.0642755419845966\t0.12140989191341\t-0.0357453508506416\t-0.0226383066061357\t0.0368248043426476\t-0.0528917387855805\t-0.00531830077284106\t-0.023209485392236\t-0.0409443982548384\t-0.0602279970133411\t-0.0532030548529174\t-0.00160243109198704\t0.01863715284029\t-0.0217948104416973\t0.383633603518805\t0.01914819572774'..b'3649\t-0.0233728418150985\t0.0834227856603043\t-0.0387212689728099\t-0.00119746419317374\t-0.0467161661199405\t0.0244932159579094\t-0.00601371599428177\t-0.0123967408954616\t-0.00869885370403252\t-0.070959838831463\t0.0253915225975361\t-0.0423681155297438\t0.066740929649166\t-0.0668503614509525\t0.0269407524752345\t-0.00233698701887992\t-0.0201744379995652\t-0.0163105807972654\t-0.0188513528988461\t0.0516066105249039\t0.0255643695407259\t-0.00547247376816918\t0.0380095529463691\t-0.00161922449024369\t0.00858485163699162\t0.0292988963958579\t0.104982968311972\t-0.0220692781820458\t-0.00891347682848018\t0.0298701808856713\t0.0240580701186076\t0.0692500526101762\t-0.0235579021404434\t0.012692833673923\t0.024384481745959\t-0.0445686540359992\t-0.0344309434556812\t-0.0653036296999577\t-0.00756293282010798\t0.0138035280951421\t-0.0401364818618468\t0.0104302840705066\t-0.0029088468277865\t-0.0107767643122829\t-0.0235683437736187\t-0.0187729075754224\t-0.00438232671397056\t0.0274281274458124\t-0.0161699796261494\t-0.0663784727396207\t-0.0418881391559233\t-0.00481668371863172\t-0.0410566291307913\t-0.0338933283159119\t-0.00937778544759433\t-0.00130300323801841\t0.055052404069415\t0.0272989236756222\t0.0492511287519239\t-0.0165008213350996\t-0.0312008344819986\t-0.00753191175962406\t-0.0396153451744063\t-0.0298181399621084\t-0.0229440562598277\t-0.0169497052694034\t-0.0325792313249729\t-0.0380530713878139\t-0.0411018357617907\t-0.0254266599012937\t-0.0260672895278166\t-0.0468103386412334\t-0.0106831380355216\t-0.0355239985602283\t0.0670808845954332\t0.0864890488487473\t0.0182130205667072\t0.0195555690436822\t-0.0188543557044614\t-0.0218374355302735\t-0.0366043997351607\t-0.00448192595311599\t0.076643616882784\t-0.0232722746254099\t-0.0167757161315926\t-0.0442757686729965\t-0.0145105434084647\t-0.0231377142871792\t0.0194959356188783\t0.0265972315541988\t-0.00750674556397457\t-0.0371290973112671\t0.0503969998560216\t-0.0141354039748982\t-0.000521853016547704\t0.0181551391924321\t0.0214674809714153\t-0.0284042369259341\t0.0460454341209673\t0.0206236337276342\t-0.0188020902890624\t-0.0587590909834051\t-0.0393496851983061\t0.00437007263095457\t0.0645800193407657\t-0.0269265729971044\t5.9070976589393e-05\t0.0504497510397723\t0.0438085910871115\t0.0278799044932392\t-0.00863443044386983\t0.0755656081415106\t0.0106476728559812\t-0.0296558738996254\t-0.00733279597293942\t0.01006873770634\t-0.0220541759287482\t-0.0141782442498728\t-0.0233657696102539\t-0.0328518687984239\t-0.00936872856338004\t0.0149850174005931\t-0.035493947797972\t-0.036104975662423\t-0.0415617421013975\t0.0318899292452593\t0.0424448214848687\t-0.00321264017144473\t-0.0110393242605037\t0.0183395650996833\t0.0761426953071529\t-0.022676664326313\t-0.0668408756679104\t-0.0424299330059493\t-0.049098679187079\t0.0106415032318233\t0.0038643454095471\t0.0087400742705149\t-0.0307173099525659\t-0.07243410214539\t0.0492952787007339\t0.0234348367260992\t0.0164755189805327\t-0.0329975012127531\t0.0694786615024704\t-0.0292641854007627\t-0.0160578776015711\t-0.0230316144888508\t0.0678973960666085\t0.0155229228376435\t-0.0217775131318073\t-0.00986460901480166\t0.0246133164419596\t-0.0277774275432409\t0.0083351721164035\t-0.0299566932514561\t-0.043083585895876\t0.0210162703188898\t-0.0144482773953615\t-0.00518019778442406\t0.0143234540764356\t0.00735928852936764\t-0.0288547232243484\t-0.0174413011665707\t0.00387800949931972\t0.0696187314048713\t-0.0280005866478152\t-0.0156130158973414\t-0.0150552391028106\t0.0186348738503162\t0.0293425545878693\t-0.0227888101434769\t0.0315757811742982\t-0.0575068044437745\t-0.0252743630110878\t-0.0272391916688631\t0.0108906973054024\t-0.0582951084454989\t-0.000913609431998681\t-0.0238016357160757\t-0.0168456632230296\t0.0900250388436311\t0.0342851981214489\t0.0225537166770209\t0.102274546570816\t-0.0423995637152744\t-0.0261467224738865\t0.0530781618319353\t-0.0108491267323352\t-0.0231280763528897\t-0.0485786625441491\t-0.0215640816919691\t-0.0490896912955414\t-0.00928793207982676\t0.0618281456247728\t-0.012038421287514\t0.0177785392095407\t-0.0385111174184926\t-0.009985449667779\t0.0048707703379256\t-0.0263785881737289\t-0.0716460773103446\t0.945922162622299\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/map.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/map.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| b'@@ -0,0 +1,50001 @@\n+SNP\tChr\tPos\r\n+SNP1\t1\t3102\r\n+SNP2\t1\t4648\r\n+SNP3\t1\t7601\r\n+SNP4\t1\t10449\r\n+SNP5\t1\t13045\r\n+SNP6\t1\t26288\r\n+SNP7\t1\t32807\r\n+SNP8\t1\t34599\r\n+SNP9\t1\t35856\r\n+SNP10\t1\t39751\r\n+SNP11\t1\t44567\r\n+SNP12\t1\t46912\r\n+SNP13\t1\t48181\r\n+SNP14\t1\t49080\r\n+SNP15\t1\t51706\r\n+SNP16\t1\t51878\r\n+SNP17\t1\t52202\r\n+SNP18\t1\t55684\r\n+SNP19\t1\t60772\r\n+SNP20\t1\t61266\r\n+SNP21\t1\t63759\r\n+SNP22\t1\t63915\r\n+SNP23\t1\t64149\r\n+SNP24\t1\t70933\r\n+SNP25\t1\t71348\r\n+SNP26\t1\t71868\r\n+SNP27\t1\t72138\r\n+SNP28\t1\t72924\r\n+SNP29\t1\t75481\r\n+SNP30\t1\t75721\r\n+SNP31\t1\t76879\r\n+SNP32\t1\t77127\r\n+SNP33\t1\t77140\r\n+SNP34\t1\t80216\r\n+SNP35\t1\t80374\r\n+SNP36\t1\t80400\r\n+SNP37\t1\t83177\r\n+SNP38\t1\t84144\r\n+SNP39\t1\t84558\r\n+SNP40\t1\t85561\r\n+SNP41\t1\t85860\r\n+SNP42\t1\t87060\r\n+SNP43\t1\t92866\r\n+SNP44\t1\t95018\r\n+SNP45\t1\t95225\r\n+SNP46\t1\t97473\r\n+SNP47\t1\t99001\r\n+SNP48\t1\t99456\r\n+SNP49\t1\t99757\r\n+SNP50\t1\t100314\r\n+SNP51\t1\t103453\r\n+SNP52\t1\t103654\r\n+SNP53\t1\t103963\r\n+SNP54\t1\t105064\r\n+SNP55\t1\t105282\r\n+SNP56\t1\t105894\r\n+SNP57\t1\t107585\r\n+SNP58\t1\t109472\r\n+SNP59\t1\t113076\r\n+SNP60\t1\t114737\r\n+SNP61\t1\t114915\r\n+SNP62\t1\t117199\r\n+SNP63\t1\t118441\r\n+SNP64\t1\t120328\r\n+SNP65\t1\t120456\r\n+SNP66\t1\t122503\r\n+SNP67\t1\t125155\r\n+SNP68\t1\t132773\r\n+SNP69\t1\t133144\r\n+SNP70\t1\t135519\r\n+SNP71\t1\t136211\r\n+SNP72\t1\t137545\r\n+SNP73\t1\t139991\r\n+SNP74\t1\t141556\r\n+SNP75\t1\t144832\r\n+SNP76\t1\t146220\r\n+SNP77\t1\t151942\r\n+SNP78\t1\t153104\r\n+SNP79\t1\t156619\r\n+SNP80\t1\t158640\r\n+SNP81\t1\t164553\r\n+SNP82\t1\t165386\r\n+SNP83\t1\t170122\r\n+SNP84\t1\t170319\r\n+SNP85\t1\t172182\r\n+SNP86\t1\t173644\r\n+SNP87\t1\t177611\r\n+SNP88\t1\t179209\r\n+SNP89\t1\t180123\r\n+SNP90\t1\t183655\r\n+SNP91\t1\t185783\r\n+SNP92\t1\t198692\r\n+SNP93\t1\t200612\r\n+SNP94\t1\t201161\r\n+SNP95\t1\t201650\r\n+SNP96\t1\t201754\r\n+SNP97\t1\t202392\r\n+SNP98\t1\t205710\r\n+SNP99\t1\t207154\r\n+SNP100\t1\t207187\r\n+SNP101\t1\t212440\r\n+SNP102\t1\t216654\r\n+SNP103\t1\t221441\r\n+SNP104\t1\t222932\r\n+SNP105\t1\t223306\r\n+SNP106\t1\t239739\r\n+SNP107\t1\t240440\r\n+SNP108\t1\t241299\r\n+SNP109\t1\t242504\r\n+SNP110\t1\t242603\r\n+SNP111\t1\t242678\r\n+SNP112\t1\t242911\r\n+SNP113\t1\t245788\r\n+SNP114\t1\t245871\r\n+SNP115\t1\t246084\r\n+SNP116\t1\t247357\r\n+SNP117\t1\t247683\r\n+SNP118\t1\t249083\r\n+SNP119\t1\t250718\r\n+SNP120\t1\t251145\r\n+SNP121\t1\t252841\r\n+SNP122\t1\t253289\r\n+SNP123\t1\t256194\r\n+SNP124\t1\t265099\r\n+SNP125\t1\t265611\r\n+SNP126\t1\t266159\r\n+SNP127\t1\t267458\r\n+SNP128\t1\t268235\r\n+SNP129\t1\t269680\r\n+SNP130\t1\t274180\r\n+SNP131\t1\t274842\r\n+SNP132\t1\t275504\r\n+SNP133\t1\t275624\r\n+SNP134\t1\t277613\r\n+SNP135\t1\t278674\r\n+SNP136\t1\t279531\r\n+SNP137\t1\t284444\r\n+SNP138\t1\t284989\r\n+SNP139\t1\t287160\r\n+SNP140\t1\t288548\r\n+SNP141\t1\t289260\r\n+SNP142\t1\t289857\r\n+SNP143\t1\t296776\r\n+SNP144\t1\t304984\r\n+SNP145\t1\t305362\r\n+SNP146\t1\t306663\r\n+SNP147\t1\t310925\r\n+SNP148\t1\t311769\r\n+SNP149\t1\t321615\r\n+SNP150\t1\t333482\r\n+SNP151\t1\t335469\r\n+SNP152\t1\t342966\r\n+SNP153\t1\t350432\r\n+SNP154\t1\t351689\r\n+SNP155\t1\t353230\r\n+SNP156\t1\t355316\r\n+SNP157\t1\t355976\r\n+SNP158\t1\t357001\r\n+SNP159\t1\t360888\r\n+SNP160\t1\t363817\r\n+SNP161\t1\t367218\r\n+SNP162\t1\t367914\r\n+SNP163\t1\t372172\r\n+SNP164\t1\t372809\r\n+SNP165\t1\t375499\r\n+SNP166\t1\t381773\r\n+SNP167\t1\t389853\r\n+SNP168\t1\t390211\r\n+SNP169\t1\t394901\r\n+SNP170\t1\t399920\r\n+SNP171\t1\t400145\r\n+SNP172\t1\t408649\r\n+SNP173\t1\t410242\r\n+SNP174\t1\t411663\r\n+SNP175\t1\t413800\r\n+SNP176\t1\t414360\r\n+SNP177\t1\t415588\r\n+SNP178\t1\t416491\r\n+SNP179\t1\t419133\r\n+SNP180\t1\t420866\r\n+SNP181\t1\t421839\r\n+SNP182\t1\t422239\r\n+SNP183\t1\t422491\r\n+SNP184\t1\t426854\r\n+SNP185\t1\t428014\r\n+SNP186\t1\t430462\r\n+SNP187\t1\t434620\r\n+SNP188\t1\t434762\r\n+SNP189\t1\t435420\r\n+SNP190\t1\t435990\r\n+SNP191\t1\t436113\r\n+SNP192\t1\t436777\r\n+SNP193\t1\t437174\r\n+SNP194\t1\t437189\r\n+SNP195\t1\t439785\r\n+SNP196\t1\t441892\r\n+SNP197\t1\t442244\r\n+SNP198\t1\t442725\r\n+SNP199\t1\t446756\r\n+SNP200\t1\t449510\r\n+SNP201\t1\t453120\r\n+SNP202\t1\t455996\r\n+SNP203\t1\t456097\r\n+SNP204\t1\t459037\r\n+SNP205\t1\t459541\r\n+SNP206\t1\t463226\r\n+SNP207\t1\t467200\r\n+SNP208\t1\t470901\r\n+SNP209\t1\t472498\r\n+SNP210\t1\t477287\r\n+SNP211\t1\t478552\r\n+SNP212\t1\t480540\r\n+SNP213\t1\t480884\r\n+SNP214\t1\t482966\r\n+SNP215\t1\t486440\r\n+SNP216\t1\t489844\r\n+SNP217\t1\t491701\r\n+SNP218\t1\t495164\r\n+SNP219\t1\t495626\r\n+SNP220\t1\t496730\r\n+SNP221\t1\t497991\r\n+SNP222\t1\t501627\r\n+SNP223\t1\t503930\r\n+SNP224\t1\t506716\r\n+SNP225\t1\t514679\r\n+SNP226\t1\t519069\r\n+SNP227\t1\t520997\r\n+SNP228\t1\t522600\r\n+SNP229\t1\t525870\r\n+SN'..b'9819\t5\t26549125\r\n+SNP49820\t5\t26549597\r\n+SNP49821\t5\t26551006\r\n+SNP49822\t5\t26554882\r\n+SNP49823\t5\t26555983\r\n+SNP49824\t5\t26556155\r\n+SNP49825\t5\t26558038\r\n+SNP49826\t5\t26559749\r\n+SNP49827\t5\t26561023\r\n+SNP49828\t5\t26562024\r\n+SNP49829\t5\t26565918\r\n+SNP49830\t5\t26570410\r\n+SNP49831\t5\t26571213\r\n+SNP49832\t5\t26580132\r\n+SNP49833\t5\t26581355\r\n+SNP49834\t5\t26584850\r\n+SNP49835\t5\t26585296\r\n+SNP49836\t5\t26588498\r\n+SNP49837\t5\t26589208\r\n+SNP49838\t5\t26589890\r\n+SNP49839\t5\t26590029\r\n+SNP49840\t5\t26590223\r\n+SNP49841\t5\t26593944\r\n+SNP49842\t5\t26596104\r\n+SNP49843\t5\t26600249\r\n+SNP49844\t5\t26602027\r\n+SNP49845\t5\t26604232\r\n+SNP49846\t5\t26609233\r\n+SNP49847\t5\t26609908\r\n+SNP49848\t5\t26611182\r\n+SNP49849\t5\t26611232\r\n+SNP49850\t5\t26612624\r\n+SNP49851\t5\t26612661\r\n+SNP49852\t5\t26612903\r\n+SNP49853\t5\t26615384\r\n+SNP49854\t5\t26615693\r\n+SNP49855\t5\t26616176\r\n+SNP49856\t5\t26617672\r\n+SNP49857\t5\t26618891\r\n+SNP49858\t5\t26621940\r\n+SNP49859\t5\t26624572\r\n+SNP49860\t5\t26630473\r\n+SNP49861\t5\t26630727\r\n+SNP49862\t5\t26631361\r\n+SNP49863\t5\t26631724\r\n+SNP49864\t5\t26632924\r\n+SNP49865\t5\t26633029\r\n+SNP49866\t5\t26635213\r\n+SNP49867\t5\t26635287\r\n+SNP49868\t5\t26637323\r\n+SNP49869\t5\t26637775\r\n+SNP49870\t5\t26641664\r\n+SNP49871\t5\t26641682\r\n+SNP49872\t5\t26644781\r\n+SNP49873\t5\t26645924\r\n+SNP49874\t5\t26646322\r\n+SNP49875\t5\t26646623\r\n+SNP49876\t5\t26648090\r\n+SNP49877\t5\t26649074\r\n+SNP49878\t5\t26650159\r\n+SNP49879\t5\t26653097\r\n+SNP49880\t5\t26653398\r\n+SNP49881\t5\t26656123\r\n+SNP49882\t5\t26656310\r\n+SNP49883\t5\t26657906\r\n+SNP49884\t5\t26663663\r\n+SNP49885\t5\t26668988\r\n+SNP49886\t5\t26672526\r\n+SNP49887\t5\t26676933\r\n+SNP49888\t5\t26678074\r\n+SNP49889\t5\t26679937\r\n+SNP49890\t5\t26680632\r\n+SNP49891\t5\t26680762\r\n+SNP49892\t5\t26693492\r\n+SNP49893\t5\t26693663\r\n+SNP49894\t5\t26694625\r\n+SNP49895\t5\t26695258\r\n+SNP49896\t5\t26697632\r\n+SNP49897\t5\t26698713\r\n+SNP49898\t5\t26699300\r\n+SNP49899\t5\t26701797\r\n+SNP49900\t5\t26712065\r\n+SNP49901\t5\t26712447\r\n+SNP49902\t5\t26714000\r\n+SNP49903\t5\t26719773\r\n+SNP49904\t5\t26720488\r\n+SNP49905\t5\t26722441\r\n+SNP49906\t5\t26722955\r\n+SNP49907\t5\t26730223\r\n+SNP49908\t5\t26730417\r\n+SNP49909\t5\t26730580\r\n+SNP49910\t5\t26730932\r\n+SNP49911\t5\t26737771\r\n+SNP49912\t5\t26741935\r\n+SNP49913\t5\t26744423\r\n+SNP49914\t5\t26750630\r\n+SNP49915\t5\t26751061\r\n+SNP49916\t5\t26753820\r\n+SNP49917\t5\t26753974\r\n+SNP49918\t5\t26754400\r\n+SNP49919\t5\t26755093\r\n+SNP49920\t5\t26757031\r\n+SNP49921\t5\t26757531\r\n+SNP49922\t5\t26757861\r\n+SNP49923\t5\t26761492\r\n+SNP49924\t5\t26764320\r\n+SNP49925\t5\t26768820\r\n+SNP49926\t5\t26770148\r\n+SNP49927\t5\t26772713\r\n+SNP49928\t5\t26775126\r\n+SNP49929\t5\t26775142\r\n+SNP49930\t5\t26781770\r\n+SNP49931\t5\t26785376\r\n+SNP49932\t5\t26801288\r\n+SNP49933\t5\t26806802\r\n+SNP49934\t5\t26813488\r\n+SNP49935\t5\t26814639\r\n+SNP49936\t5\t26817151\r\n+SNP49937\t5\t26827902\r\n+SNP49938\t5\t26832589\r\n+SNP49939\t5\t26833631\r\n+SNP49940\t5\t26836247\r\n+SNP49941\t5\t26837168\r\n+SNP49942\t5\t26845045\r\n+SNP49943\t5\t26845380\r\n+SNP49944\t5\t26846791\r\n+SNP49945\t5\t26846997\r\n+SNP49946\t5\t26848326\r\n+SNP49947\t5\t26848479\r\n+SNP49948\t5\t26848902\r\n+SNP49949\t5\t26853769\r\n+SNP49950\t5\t26854942\r\n+SNP49951\t5\t26856026\r\n+SNP49952\t5\t26865256\r\n+SNP49953\t5\t26866572\r\n+SNP49954\t5\t26867935\r\n+SNP49955\t5\t26871770\r\n+SNP49956\t5\t26872990\r\n+SNP49957\t5\t26873413\r\n+SNP49958\t5\t26876546\r\n+SNP49959\t5\t26877047\r\n+SNP49960\t5\t26877535\r\n+SNP49961\t5\t26880272\r\n+SNP49962\t5\t26884403\r\n+SNP49963\t5\t26885612\r\n+SNP49964\t5\t26886314\r\n+SNP49965\t5\t26888824\r\n+SNP49966\t5\t26890822\r\n+SNP49967\t5\t26898401\r\n+SNP49968\t5\t26898888\r\n+SNP49969\t5\t26900331\r\n+SNP49970\t5\t26900828\r\n+SNP49971\t5\t26912662\r\n+SNP49972\t5\t26913253\r\n+SNP49973\t5\t26914777\r\n+SNP49974\t5\t26915718\r\n+SNP49975\t5\t26917075\r\n+SNP49976\t5\t26918397\r\n+SNP49977\t5\t26927110\r\n+SNP49978\t5\t26928159\r\n+SNP49979\t5\t26928701\r\n+SNP49980\t5\t26930248\r\n+SNP49981\t5\t26930710\r\n+SNP49982\t5\t26931589\r\n+SNP49983\t5\t26931709\r\n+SNP49984\t5\t26932178\r\n+SNP49985\t5\t26935454\r\n+SNP49986\t5\t26936463\r\n+SNP49987\t5\t26937752\r\n+SNP49988\t5\t26937854\r\n+SNP49989\t5\t26939596\r\n+SNP49990\t5\t26940681\r\n+SNP49991\t5\t26941036\r\n+SNP49992\t5\t26941289\r\n+SNP49993\t5\t26941789\r\n+SNP49994\t5\t26947445\r\n+SNP49995\t5\t26960324\r\n+SNP49996\t5\t26964012\r\n+SNP49997\t5\t26964614\r\n+SNP49998\t5\t26966149\r\n+SNP49999\t5\t26971332\r\n+SNP50000\t5\t26973018\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/output.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| b'@@ -0,0 +1,50001 @@\n+Trait\tMarker_name\tChr\tPos\tPvalue\n+traitname\tSNP1\t1\t3102\t0.469258063992063\n+traitname\tSNP2\t1\t4648\t0.101620713301417\n+traitname\tSNP3\t1\t7601\t0.730848388413714\n+traitname\tSNP4\t1\t10449\t0.782552899983113\n+traitname\tSNP5\t1\t13045\t0.0219306740805478\n+traitname\tSNP6\t1\t26288\t0.88617546778878\n+traitname\tSNP7\t1\t32807\t0.972040107240628\n+traitname\tSNP8\t1\t34599\t0.498651882283497\n+traitname\tSNP9\t1\t35856\t0.111605610470404\n+traitname\tSNP10\t1\t39751\t0.423226429567156\n+traitname\tSNP11\t1\t44567\t0.481335913032479\n+traitname\tSNP12\t1\t46912\t0.108124737530005\n+traitname\tSNP13\t1\t48181\t0.695641883380419\n+traitname\tSNP14\t1\t49080\t0.856375273530245\n+traitname\tSNP15\t1\t51706\t0.968487529146441\n+traitname\tSNP16\t1\t51878\t0.693301520986122\n+traitname\tSNP17\t1\t52202\t0.49522970348474\n+traitname\tSNP18\t1\t55684\t0.75059479811401\n+traitname\tSNP19\t1\t60772\t0.518806774519514\n+traitname\tSNP20\t1\t61266\t0.931984181222639\n+traitname\tSNP21\t1\t63759\t0.643701861734718\n+traitname\tSNP22\t1\t63915\t0.272594443915554\n+traitname\tSNP23\t1\t64149\t0.663966971107014\n+traitname\tSNP24\t1\t70933\t0.90694100361626\n+traitname\tSNP25\t1\t71348\t0.154206636081275\n+traitname\tSNP26\t1\t71868\t0.446613665755518\n+traitname\tSNP27\t1\t72138\t0.491693690832368\n+traitname\tSNP28\t1\t72924\t0.346078261207293\n+traitname\tSNP29\t1\t75481\t0.131428010841624\n+traitname\tSNP30\t1\t75721\t0.751727481267894\n+traitname\tSNP31\t1\t76879\t0.261144989676999\n+traitname\tSNP32\t1\t77127\t0.449192097599046\n+traitname\tSNP33\t1\t77140\t0.506568467056256\n+traitname\tSNP34\t1\t80216\t0.366091908315112\n+traitname\tSNP35\t1\t80374\t0.350183343316253\n+traitname\tSNP36\t1\t80400\t0.614728775987361\n+traitname\tSNP37\t1\t83177\t0.207711200609552\n+traitname\tSNP38\t1\t84144\t0.683185043910651\n+traitname\tSNP39\t1\t84558\t0.753550329989218\n+traitname\tSNP40\t1\t85561\t0.76889645811777\n+traitname\tSNP41\t1\t85860\t0.288849105949473\n+traitname\tSNP42\t1\t87060\t0.598751140714578\n+traitname\tSNP43\t1\t92866\t0.861787634523111\n+traitname\tSNP44\t1\t95018\t0.269429297048646\n+traitname\tSNP45\t1\t95225\t0.772581810601898\n+traitname\tSNP46\t1\t97473\t0.648431702280633\n+traitname\tSNP47\t1\t99001\t0.600529159325437\n+traitname\tSNP48\t1\t99456\t0.709504544599765\n+traitname\tSNP49\t1\t99757\t0.339299730138436\n+traitname\tSNP50\t1\t100314\t0.77122977426066\n+traitname\tSNP51\t1\t103453\t0.320153193461896\n+traitname\tSNP52\t1\t103654\t0.100871264553898\n+traitname\tSNP53\t1\t103963\t0.767323169723174\n+traitname\tSNP54\t1\t105064\t0.00715808977364939\n+traitname\tSNP55\t1\t105282\t0.250266944100426\n+traitname\tSNP56\t1\t105894\t0.102021249052409\n+traitname\tSNP57\t1\t107585\t0.96622032714115\n+traitname\tSNP58\t1\t109472\t0.213782862586124\n+traitname\tSNP59\t1\t113076\t0.143226320707381\n+traitname\tSNP60\t1\t114737\t0.924013276218658\n+traitname\tSNP61\t1\t114915\t0.745546324942985\n+traitname\tSNP62\t1\t117199\t0.645888811224891\n+traitname\tSNP63\t1\t118441\t0.495512048847985\n+traitname\tSNP64\t1\t120328\t0.785061371868697\n+traitname\tSNP65\t1\t120456\t0.700632087371393\n+traitname\tSNP66\t1\t122503\t0.169284280137093\n+traitname\tSNP67\t1\t125155\t0.217071312332108\n+traitname\tSNP68\t1\t132773\t0.159886400029993\n+traitname\tSNP69\t1\t133144\t0.831698817765358\n+traitname\tSNP70\t1\t135519\t0.249360550510533\n+traitname\tSNP71\t1\t136211\t0.518356218314296\n+traitname\tSNP72\t1\t137545\t0.299154333934361\n+traitname\tSNP73\t1\t139991\t0.947058371085541\n+traitname\tSNP74\t1\t141556\t0.756108502265097\n+traitname\tSNP75\t1\t144832\t0.271154291228045\n+traitname\tSNP76\t1\t146220\t0.88236837881742\n+traitname\tSNP77\t1\t151942\t0.170812108410772\n+traitname\tSNP78\t1\t153104\t0.365914352292149\n+traitname\tSNP79\t1\t156619\t0.767187083780078\n+traitname\tSNP80\t1\t158640\t0.860292030226417\n+traitname\tSNP81\t1\t164553\t0.167719947407506\n+traitname\tSNP82\t1\t165386\t0.330794226336004\n+traitname\tSNP83\t1\t170122\t0.692452416302984\n+traitname\tSNP84\t1\t170319\t0.0803059486518452\n+traitname\tSNP85\t1\t172182\t0.582547978499917\n+traitname\tSNP86\t1\t173644\t0.664355087833467\n+traitname\tSNP87\t1\t177611\t0.574903641394203\n+traitname\tSNP88\t1\t179209\t0.417450551553441\n+traitname\tSNP89\t1\t180123\t0.999391968979726\n+traitname\tSNP90\t1\t183655\t0.686904731263637\n+traitname\tSNP91\t1\t185783\t0.30768721041746\n+traitn'..b'\t26755093\t0.369861981625335\n+traitname\tSNP49920\t5\t26757031\t0.459339311458927\n+traitname\tSNP49921\t5\t26757531\t0.0909352373147266\n+traitname\tSNP49922\t5\t26757861\t0.589237922890095\n+traitname\tSNP49923\t5\t26761492\t0.968819604168695\n+traitname\tSNP49924\t5\t26764320\t0.186929594593407\n+traitname\tSNP49925\t5\t26768820\t0.708434665277452\n+traitname\tSNP49926\t5\t26770148\t0.473799947083281\n+traitname\tSNP49927\t5\t26772713\t0.724392066172182\n+traitname\tSNP49928\t5\t26775126\t0.661226113373206\n+traitname\tSNP49929\t5\t26775142\t0.181303027584578\n+traitname\tSNP49930\t5\t26781770\t0.235863227652609\n+traitname\tSNP49931\t5\t26785376\t0.717018050219037\n+traitname\tSNP49932\t5\t26801288\t0.703597011068272\n+traitname\tSNP49933\t5\t26806802\t0.484710342544411\n+traitname\tSNP49934\t5\t26813488\t0.0479663598353848\n+traitname\tSNP49935\t5\t26814639\t0.596910734540523\n+traitname\tSNP49936\t5\t26817151\t0.0312083095983475\n+traitname\tSNP49937\t5\t26827902\t0.59525405631134\n+traitname\tSNP49938\t5\t26832589\t0.39815533460193\n+traitname\tSNP49939\t5\t26833631\t0.891321223937739\n+traitname\tSNP49940\t5\t26836247\t0.895914754574338\n+traitname\tSNP49941\t5\t26837168\t0.083554995408947\n+traitname\tSNP49942\t5\t26845045\t0.676045462685435\n+traitname\tSNP49943\t5\t26845380\t0.153356073373809\n+traitname\tSNP49944\t5\t26846791\t0.313657704645926\n+traitname\tSNP49945\t5\t26846997\t0.345703221263434\n+traitname\tSNP49946\t5\t26848326\t0.301080087365683\n+traitname\tSNP49947\t5\t26848479\t0.162055136624132\n+traitname\tSNP49948\t5\t26848902\t0.900402672935511\n+traitname\tSNP49949\t5\t26853769\t0.0429262765494099\n+traitname\tSNP49950\t5\t26854942\t0.532673542562009\n+traitname\tSNP49951\t5\t26856026\t0.615437137895911\n+traitname\tSNP49952\t5\t26865256\t0.175912444752761\n+traitname\tSNP49953\t5\t26866572\t0.742170597671728\n+traitname\tSNP49954\t5\t26867935\t0.288836612497629\n+traitname\tSNP49955\t5\t26871770\t0.98613663441062\n+traitname\tSNP49956\t5\t26872990\t0.00440990480672004\n+traitname\tSNP49957\t5\t26873413\t0.487722875010919\n+traitname\tSNP49958\t5\t26876546\t0.0270410359639908\n+traitname\tSNP49959\t5\t26877047\t0.10594425326491\n+traitname\tSNP49960\t5\t26877535\t0.583857740091527\n+traitname\tSNP49961\t5\t26880272\t0.0525001910950511\n+traitname\tSNP49962\t5\t26884403\t0.499654023075892\n+traitname\tSNP49963\t5\t26885612\t0.0154073869214085\n+traitname\tSNP49964\t5\t26886314\t0.383457389725903\n+traitname\tSNP49965\t5\t26888824\t0.363597013001999\n+traitname\tSNP49966\t5\t26890822\t0.563241800147088\n+traitname\tSNP49967\t5\t26898401\t0.64076216246539\n+traitname\tSNP49968\t5\t26898888\t0.165271494517212\n+traitname\tSNP49969\t5\t26900331\t0.922340602595221\n+traitname\tSNP49970\t5\t26900828\t0.170326751757837\n+traitname\tSNP49971\t5\t26912662\t0.888946047997379\n+traitname\tSNP49972\t5\t26913253\t0.871955104159578\n+traitname\tSNP49973\t5\t26914777\t0.424041808890503\n+traitname\tSNP49974\t5\t26915718\t0.823006879593037\n+traitname\tSNP49975\t5\t26917075\t0.375798520392157\n+traitname\tSNP49976\t5\t26918397\t0.213018327939401\n+traitname\tSNP49977\t5\t26927110\t0.97026274750341\n+traitname\tSNP49978\t5\t26928159\t0.786528008378223\n+traitname\tSNP49979\t5\t26928701\t0.55637341896506\n+traitname\tSNP49980\t5\t26930248\t0.928984251562843\n+traitname\tSNP49981\t5\t26930710\t0.13986357946646\n+traitname\tSNP49982\t5\t26931589\t0.300800695308537\n+traitname\tSNP49983\t5\t26931709\t0.430249428092055\n+traitname\tSNP49984\t5\t26932178\t0.989977997472187\n+traitname\tSNP49985\t5\t26935454\t0.00692951153926428\n+traitname\tSNP49986\t5\t26936463\t0.980332090447985\n+traitname\tSNP49987\t5\t26937752\t0.858543442379832\n+traitname\tSNP49988\t5\t26937854\t0.801454701969988\n+traitname\tSNP49989\t5\t26939596\t0.171316978188681\n+traitname\tSNP49990\t5\t26940681\t0.145664750742478\n+traitname\tSNP49991\t5\t26941036\t0.447315064955819\n+traitname\tSNP49992\t5\t26941289\t0.155508115777885\n+traitname\tSNP49993\t5\t26941789\t0.167886917739672\n+traitname\tSNP49994\t5\t26947445\t0.787787410972909\n+traitname\tSNP49995\t5\t26960324\t0.20210766991913\n+traitname\tSNP49996\t5\t26964012\t0.815753119722407\n+traitname\tSNP49997\t5\t26964614\t0.291933144356384\n+traitname\tSNP49998\t5\t26966149\t0.526026466113253\n+traitname\tSNP49999\t5\t26971332\t0.0125514745273496\n+traitname\tSNP50000\t5\t26973018\t0.49073221323915\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/phenot.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/phenot.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| b'@@ -0,0 +1,501 @@\n+Ind_id\tPhenot1\tPhenot2\r\n+Ind1\t-2.9985936006411\t-2.9985936006411\r\n+Ind2\t-2.68669426456267\t-2.68669426456267\r\n+Ind3\t-18.2181577678362\t-18.2181577678362\r\n+Ind4\t-12.5401016500159\t-12.5401016500159\r\n+Ind5\t7.72670043773802\t7.72670043773802\r\n+Ind6\t5.25010246801007\t5.25010246801007\r\n+Ind7\t0.417402126827894\t0.417402126827894\r\n+Ind8\t10.6214497655125\t10.6214497655125\r\n+Ind9\t-18.3014674338768\t-18.3014674338768\r\n+Ind10\t-18.6098675031914\tNA\r\n+Ind11\t-4.03477186214533\t-4.03477186214533\r\n+Ind12\t-15.9882829569252\t-15.9882829569252\r\n+Ind13\t-1.04632067604619\t-1.04632067604619\r\n+Ind14\t-5.66838909636386\t-5.66838909636386\r\n+Ind15\t2.58509745805752\t2.58509745805752\r\n+Ind16\t-3.24187957417854\t-3.24187957417854\r\n+Ind17\t-23.2729677083008\t-23.2729677083008\r\n+Ind18\t6.96435357361898\t6.96435357361898\r\n+Ind19\t-8.86432852665433\t-8.86432852665433\r\n+Ind20\t-2.43021151820351\t-2.43021151820351\r\n+Ind21\t-11.4044395258052\t-11.4044395258052\r\n+Ind22\t4.77444855571111\t4.77444855571111\r\n+Ind23\t18.2425783461208\tNA\r\n+Ind24\t-16.1003230103367\t-16.1003230103367\r\n+Ind25\t3.48578405180521\t3.48578405180521\r\n+Ind26\t-13.6317796013066\t-13.6317796013066\r\n+Ind27\t-17.8299036740888\t-17.8299036740888\r\n+Ind28\t12.8760218262878\t12.8760218262878\r\n+Ind29\t0.982337452944073\t0.982337452944073\r\n+Ind30\t6.89441228302578\tNA\r\n+Ind31\t-9.79070987759412\t-9.79070987759412\r\n+Ind32\t-7.33335887602666\t-7.33335887602666\r\n+Ind33\t-9.38112263525387\t-9.38112263525387\r\n+Ind34\t-17.418514681802\t-17.418514681802\r\n+Ind35\t-1.36926657800065\tNA\r\n+Ind36\t-15.8925185459182\t-15.8925185459182\r\n+Ind37\t-16.6664415036233\t-16.6664415036233\r\n+Ind38\t-17.6113176227666\t-17.6113176227666\r\n+Ind39\t8.09960655518536\t8.09960655518536\r\n+Ind40\t-14.9600772642048\t-14.9600772642048\r\n+Ind41\t-2.46258520207899\t-2.46258520207899\r\n+Ind42\t3.52692956060259\t3.52692956060259\r\n+Ind43\t-2.56849421725468\tNA\r\n+Ind44\t0.462204831714866\t0.462204831714866\r\n+Ind45\t-3.38429200532215\t-3.38429200532215\r\n+Ind46\t-9.02868604003372\t-9.02868604003372\r\n+Ind47\t5.19312365030726\t5.19312365030726\r\n+Ind48\t11.362118479605\t11.362118479605\r\n+Ind49\t-4.7016396675323\t-4.7016396675323\r\n+Ind50\t-6.49221835766938\t-6.49221835766938\r\n+Ind51\t8.74611923535598\t8.74611923535598\r\n+Ind52\t20.0594288752544\t20.0594288752544\r\n+Ind53\t-1.11475618282811\t-1.11475618282811\r\n+Ind54\t-10.6088037144513\t-10.6088037144513\r\n+Ind55\t-1.02045642033812\t-1.02045642033812\r\n+Ind56\t-9.96139275677573\t-9.96139275677573\r\n+Ind57\t3.07304203895751\t3.07304203895751\r\n+Ind58\t-12.1924312705288\t-12.1924312705288\r\n+Ind59\t3.05865100310952\t3.05865100310952\r\n+Ind60\t7.09555942891504\t7.09555942891504\r\n+Ind61\t-11.8186149500922\t-11.8186149500922\r\n+Ind62\t-1.65105172480334\t-1.65105172480334\r\n+Ind63\t1.30018923413375\t1.30018923413375\r\n+Ind64\t-13.4333408135065\t-13.4333408135065\r\n+Ind65\t-7.02813352566701\t-7.02813352566701\r\n+Ind66\t13.8480146064731\t13.8480146064731\r\n+Ind67\t-0.984981720446856\t-0.984981720446856\r\n+Ind68\t-16.5795212756351\t-16.5795212756351\r\n+Ind69\t4.17245531503865\t4.17245531503865\r\n+Ind70\t-6.42799646672008\t-6.42799646672008\r\n+Ind71\t13.5192944193829\t13.5192944193829\r\n+Ind72\t-2.80402277594647\t-2.80402277594647\r\n+Ind73\t-15.0486960942153\t-15.0486960942153\r\n+Ind74\t0.0372137470215899\t0.0372137470215899\r\n+Ind75\t-0.854371678661337\t-0.854371678661337\r\n+Ind76\t11.5856619179338\t11.5856619179338\r\n+Ind77\t-26.4907768158664\t-26.4907768158664\r\n+Ind78\t-3.12931669517017\t-3.12931669517017\r\n+Ind79\t-13.760048976086\t-13.760048976086\r\n+Ind80\t-1.03667332000764\t-1.03667332000764\r\n+Ind81\t1.98118654229534\t1.98118654229534\r\n+Ind82\t-9.63094615407686\t-9.63094615407686\r\n+Ind83\t-11.6395266623602\t-11.6395266623602\r\n+Ind84\t-7.08357127398883\t-7.08357127398883\r\n+Ind85\t12.9387551173972\t12.9387551173972\r\n+Ind86\t-1.30834381449001\t-1.30834381449001\r\n+Ind87\t13.6609913266554\t13.6609913266554\r\n+Ind88\t6.39847190987683\t6.39847190987683\r\n+Ind89\t-1.34730952287177\t-1.34730952287177\r\n+Ind90\t-16.6048003152535\t-16.6048003152535\r\n+Ind91\t16.6084424113242\t16.6084424113242\r\n+Ind92\t6.64141034078026\tNA\r\n+Ind93\t-9.40744181889942\t-9.40744181889942\r\n+Ind94\t14.4360126037192\t'..b'9523371383942\t13.9523371383942\r\n+Ind410\t-21.7489986447953\t-21.7489986447953\r\n+Ind411\t3.63087289660466\t3.63087289660466\r\n+Ind412\t-7.28557138176745\t-7.28557138176745\r\n+Ind413\t9.90771261975745\t9.90771261975745\r\n+Ind414\t-6.63360770716923\t-6.63360770716923\r\n+Ind415\t-11.6462997607846\t-11.6462997607846\r\n+Ind416\t-5.16413899142312\t-5.16413899142312\r\n+Ind417\t4.18820760770602\t4.18820760770602\r\n+Ind418\t-7.99453947508391\t-7.99453947508391\r\n+Ind419\t-7.36353065076308\t-7.36353065076308\r\n+Ind420\t1.55602020780775\t1.55602020780775\r\n+Ind421\t-9.49723313250572\t-9.49723313250572\r\n+Ind422\t-6.88614454868852\t-6.88614454868852\r\n+Ind423\t-12.3248092127303\t-12.3248092127303\r\n+Ind424\t-18.671461132403\t-18.671461132403\r\n+Ind425\t1.36473273071079\tNA\r\n+Ind426\t-5.8867712610241\t-5.8867712610241\r\n+Ind427\t-9.50273004869248\t-9.50273004869248\r\n+Ind428\t-4.81943622063048\t-4.81943622063048\r\n+Ind429\t-26.2142525264157\t-26.2142525264157\r\n+Ind430\t14.8425350408079\t14.8425350408079\r\n+Ind431\t-9.16761728314715\t-9.16761728314715\r\n+Ind432\t-1.00320200216495\t-1.00320200216495\r\n+Ind433\t13.774402701527\t13.774402701527\r\n+Ind434\t-13.1693917708243\t-13.1693917708243\r\n+Ind435\t-15.6696471881704\t-15.6696471881704\r\n+Ind436\t-3.67159990332953\t-3.67159990332953\r\n+Ind437\t13.5297254120583\t13.5297254120583\r\n+Ind438\t-0.112602110411935\t-0.112602110411935\r\n+Ind439\t-13.8888752459368\t-13.8888752459368\r\n+Ind440\t-5.06357674519975\t-5.06357674519975\r\n+Ind441\t3.70142882013138\t3.70142882013138\r\n+Ind442\t-8.99091178928531\t-8.99091178928531\r\n+Ind443\t7.17805994018753\t7.17805994018753\r\n+Ind444\t-0.266253301970805\t-0.266253301970805\r\n+Ind445\t8.94133832711641\t8.94133832711641\r\n+Ind446\t11.8626382537432\t11.8626382537432\r\n+Ind447\t1.46579823769644\t1.46579823769644\r\n+Ind448\t4.97944326839486\t4.97944326839486\r\n+Ind449\t-6.12599552963873\t-6.12599552963873\r\n+Ind450\t2.02124592400588\t2.02124592400588\r\n+Ind451\t5.99343387340862\t5.99343387340862\r\n+Ind452\t5.91859404910017\t5.91859404910017\r\n+Ind453\t5.26076650830962\tNA\r\n+Ind454\t21.079384146455\t21.079384146455\r\n+Ind455\t-19.8109434066481\t-19.8109434066481\r\n+Ind456\t-1.95890830148838\t-1.95890830148838\r\n+Ind457\t-21.7572487645643\t-21.7572487645643\r\n+Ind458\t6.68427079420261\t6.68427079420261\r\n+Ind459\t-0.00848186185855004\t-0.00848186185855004\r\n+Ind460\t0.173826840582155\t0.173826840582155\r\n+Ind461\t2.91595224945219\t2.91595224945219\r\n+Ind462\t-9.29007019936385\t-9.29007019936385\r\n+Ind463\t0.381452713728563\t0.381452713728563\r\n+Ind464\t-3.89608517962285\t-3.89608517962285\r\n+Ind465\t13.158651656358\t13.158651656358\r\n+Ind466\t-2.4007151168472\t-2.4007151168472\r\n+Ind467\t13.4474547090651\t13.4474547090651\r\n+Ind468\t12.856692609112\t12.856692609112\r\n+Ind469\t-3.25976986805755\t-3.25976986805755\r\n+Ind470\t-6.35949771507445\t-6.35949771507445\r\n+Ind471\t20.3208227858756\tNA\r\n+Ind472\t5.67442559054238\t5.67442559054238\r\n+Ind473\t-4.47561172897561\t-4.47561172897561\r\n+Ind474\t-9.14418324318106\t-9.14418324318106\r\n+Ind475\t1.1887403231858\t1.1887403231858\r\n+Ind476\t-13.5862130319871\t-13.5862130319871\r\n+Ind477\t16.6879999050854\t16.6879999050854\r\n+Ind478\t-17.1389687195975\t-17.1389687195975\r\n+Ind479\t-13.551122507557\t-13.551122507557\r\n+Ind480\t-8.94167774042331\t-8.94167774042331\r\n+Ind481\t7.29409149872752\t7.29409149872752\r\n+Ind482\t-1.77917153389365\t-1.77917153389365\r\n+Ind483\t-7.09889887976419\t-7.09889887976419\r\n+Ind484\t-24.1386330395452\t-24.1386330395452\r\n+Ind485\t0.834467682805625\t0.834467682805625\r\n+Ind486\t-2.57299735289436\t-2.57299735289436\r\n+Ind487\t-21.3142343190533\t-21.3142343190533\r\n+Ind488\t-15.3242121153467\t-15.3242121153467\r\n+Ind489\t-1.26949587391629\t-1.26949587391629\r\n+Ind490\t-17.2632576427363\t-17.2632576427363\r\n+Ind491\t2.62910436373452\t2.62910436373452\r\n+Ind492\t-11.4255273107226\t-11.4255273107226\r\n+Ind493\t-8.51899200969042\t-8.51899200969042\r\n+Ind494\t-1.07660369147835\t-1.07660369147835\r\n+Ind495\t-13.876458194635\tNA\r\n+Ind496\t-16.303864296789\t-16.303864296789\r\n+Ind497\t0.539683578070243\t0.539683578070243\r\n+Ind498\t-0.51004693545679\t-0.51004693545679\r\n+Ind499\t-6.12081197198746\t-6.12081197198746\r\n+Ind500\t-9.94515897141559\t-9.94515897141559\r\n' |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/rss.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/rss.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| @@ -0,0 +1,21 @@ +expl_RSS h2_RSS unexpl_RSS +0 0.496182930458548 1 +0.0791505627350105 0.527910417939619 1 +0.151558378268239 0.525630045437601 1 +0.214131827973417 0.499760433439174 1 +0.280191425218391 0.490987854639079 1 +0.301542565504732 0.502911037355004 1 +0.335619292956412 0.501910798727557 1 +0.345227580376937 0.506758620219731 1 +0.385926966136037 0.496068425819074 1 +0.410642736432562 0.490945056417703 1 +0.410642736432562 0.490945056417703 1 +0.395262682359115 0.46240936659192 1 +0.376569083055453 0.46130750181542 1 +0.350355079873684 0.472140113860783 1 +0.322122466126443 0.465944614037348 1 +0.280191425218391 0.490987854639079 1 +0.214131827973417 0.499760433439174 1 +0.151558378268239 0.525630045437601 1 +0.0791505627350105 0.527910417939619 1 +0 0.496182930458548 1 |
| b |
| diff -r 6b7107812931 -r 380b364980f9 test-data/step_table.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/step_table.txt Mon Apr 16 08:50:05 2018 -0400 |
| b |
| @@ -0,0 +1,21 @@ +step step_ cof ncof h2 maxpval BIC extBIC +0 fwd0 NA 0 0.496182930458548 0 3672.84585166067 3672.84585166067 +1 fwd1 +SNP18379 1 0.48733249654522 6.90994712835942e-08 3649.84971304171 3671.48926961053 +2 fwd2 +SNP18233 2 0.440892640799306 6.2939290729918e-07 3631.37115359915 3673.26393237528 +3 fwd3 +SNP31696 3 0.363456131235323 2.85989830524578e-07 3611.7359200121 3673.07095077811 +4 fwd4 +SNP45571 4 0.292850678369099 1.48126081724671e-07 3590.56770904145 3670.76958765043 +5 fwd5 +SNP41396 5 0.288304572197428 3.25372556529604e-05 3579.29312332227 3677.9155226688 +6 fwd6 +SNP44453 6 0.250295506790257 1.55172102286074e-05 3566.4356003931 3683.11383736 +7 fwd7 +SNP22033 7 0.246697990021913 6.73987158924063e-05 3556.48901398701 3690.91474721022 +8 fwd8 +SNP3524 8 0.179362150117533 9.68164394045924e-05 3546.84278291086 3698.74890959993 +9 fwd9 +SNP2705 9 0.136254060056991 8.69490160099236e-05 3535.98995238008 3705.1408664577 +10 bwd0 -NA 9 0.136254060056991 8.69490160099236e-05 3535.98995238008 3705.1408664577 +11 bwd1 -SNP22033 8 0.111034464508902 0.000147354836597434 3545.32818192575 3697.23430861482 +12 bwd2 -SNP41396 7 0.135922708445825 2.09725261076646e-05 3553.87075026784 3688.29648349105 +13 bwd3 -SNP2705 6 0.187463997968951 1.11075179344164e-05 3565.98962143318 3682.66785840008 +14 bwd4 -SNP44453 5 0.212165384931804 2.53449108561542e-05 3579.37453327783 3677.99693262437 +15 bwd5 -SNP3524 4 0.292850678369099 1.48126081724671e-07 3590.56770904145 3670.76958765043 +16 bwd6 -SNP45571 3 0.363456131235323 2.85989830524578e-07 3611.7359200121 3673.07095077811 +17 bwd7 -SNP31696 2 0.440892640799306 6.2939290729918e-07 3631.37115359915 3673.26393237528 +18 bwd8 -SNP18233 1 0.48733249654522 6.90994712835942e-08 3649.84971304171 3671.48926961053 +19 bwd9 -SNP18379 0 0.496182930458548 0 3672.84585166067 3672.84585166067 |