SEQMIX

2013-03-13T17:53:04Z

Youna: /* Introduction */

SEQMIX

2013-03-13T17:49:10Z

Youna:

== Introduction ==

== Method ==

== Download ==

== Related Programs ==

SEQMIX

2013-03-13T17:48:35Z

Youna: Created page with '== Introduction == == Download == == Workflow == == Related Programs =='

== Introduction ==

== Download ==

== Workflow ==

== Related Programs ==

LibStatGen: VCF

2013-01-23T20:09:56Z

Youna:

[[Category:C++]]
[[Category:libStatGen]]
[[Category:libStatGen VCF]]

= Variant Call Format (VCF) =

The documentation on the Variant Call Format (VCF) can be found at: http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41

<span style="color:#D2691E">'''Our APIs are currently in the initial test phase, and are not yet available for public use.'''</span>

== API for Reading VCF Files ==

You will use both an <code>VcfFileReader</code> and an <code>VCfRecord</code> for reading VCF files.

=== <code>VcfFileReader</code> ===
An instance of the VcfFileReader class is used to read VCF files.

==== include file ====
<code>VcfFileReader</code> is declared in <code>VcfFileReader.h</code>, so be sure to include that file.
<source lang="cpp">
#include "VcfFileReader.h"
</source>

==== Opening the VCF File ====
<code>open</code> opens the specified file and by default throws an exception if it was not successfully opened.

There are two methods for opening a file for reading. One is standard while the other allows the specification of a subset of samples to keep.

Both methods take the filename to be read as well as a reference to the VcfHeader to be populated.

The standard method is:
<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
reader.open("vcfFileName.vcf", header);
</source>

==== Reading a Subset of Samples ====
To select only a subset of samples to keep, when opening the file also specify the name of the file containing the names of the samples to keep and the delimiter separating the sample names (default is a new line, '\n').
<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
// Subset 1 is delimited by new lines, '\n'.
reader.open("vcfFileName.vcf", header, "subsetFile1.txt", NULL, NULL);
</source>
<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
// Subset 2 is delimited by ';'
reader.open("vcfFileName.vcf", header, "subsetFile2.txt", NULL, NULL, ';');
</source>

To select a subset of samples to remove/ignore, when opening the file also specify the name of the file containing the names of the samples to remove/ignore and the delimiter separating the sample names (default is a new line, '\n').
<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
// Subset 1 is delimited by new lines, '\n'.
reader.open("vcfFileName.vcf", header, NULL, NULL, "subsetFile1.txt");
</source>
<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
// Subset 2 is delimited by ';'
reader.open("vcfFileName.vcf", header, NULL, NULL, "subsetFile2.txt", ';');
</source>

If you just have 1 sample to be excluded, you can specify it directly in the open line.
<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
// Subset 1 is delimited by new lines, '\n'.
reader.open("vcfFileName.vcf", header, NULL, "BAD_SAMPLE", NULL);
</source>

==== Reading a VCF Record ====
Once the file is open, get the next record using <code>readRecord</code>.

It returns true if a record was successfully found and false on EOF or an error. Typically on a reading error an exception is thrown.

<code>readRecord</code> takes a reference to a <code>VcfRecord</code> object as a parameter. When true is returned, the <code>VcfRecord</code> is updated with the next record.

If you only process one record at a time (are done with a record before reading the next one), you can loop until <code>readRecord</code> return false, reusing the same record for each call.
<source lang="cpp">
VcfRecord record;
while(reader.readRecord(record))
{
// Your record specific processing here.
}
// Done reading the file.
</source>

If any subsetting was specified when the file was opened, the record that is returned will have only the specified subset of samples in it.

When reading a record, you can also specify discard rules. Internal processing of <code>readRecord</code> will apply the discard rules on the record that was read and will continue reading records and not return from the method until a record is found that should not be discarded (or until the end of the file). True is returned if a record that should be kept was found, false if the end of the file was hit without finding a record to keep.

==== Specifying Discard Rules ====
===== Basic Rules =====
When specifying the discard rules, you should use the constants found at the top of VcfFileReader.h

For Example (see the file for the current set of Discard Rules and the associated values):
<source lang="cpp">
static const int DISCARD_NON_PHASED = 0x1;
static const int DISCARD_MISSING_GT = 0x2;
static const int DISCARD_FILTERED = 0x4;
static const int DISCARD_MULTIPLE_ALTS = 0x8;
</source>

You can set the Discard using <code>setDiscardRules(uint32_t)</code>. You can add additional rules to what has already been set by using <code>addDiscardRules(uint32_t)</code>.

<source lang="cpp">
reader.setDiscardRules(DISCARD_FILTERED);
reader.addDiscardRules(DISCARD_NON_PHASED);
</source>
is the same as
<source lang="cpp">
reader.setDiscardRules(DISCARD_FILTERED | DISCARD_NON_PHASED);
</source>
and discards reads that do not have <code>PASS</code> in the <code>FILTER</code> field and reads that have a genotype that is not phased or have no <code>GT</code> in the <code>FORMAT</code> fields.

===== Minimum Alternate Allele Count =====
To Discard any records without a minimum number of alternate alleles, use:
<source lang="cpp">
VcfFileReader::addDiscardMinAltAlleleCount(int32_t minAltAlleleCount, VcfSubsetSamples* subset)
</source>

The <code>minAltAlleleCount</code> parameter is the minimum number of alternate alleles found in the specified subset (if specified) in order for the record to be kept.

The <code>VcfSubsetSamples* subset</code> parameter is a pointer to the subset of samples that you want to include when counting the number of alternate alleles. If all samples that are read/kept are to be included, NULL should be passed in.

See [[#Handling a Subset of Samples|Handling a Subset of Samples]] for how to use <code>VcfSubsetSamples</code>.

Use the following method to remove the DiscardMinAltAlleleCount rule:
<source lang="cpp">
VcfFileReader::rmDiscardMinAltAlleleCount()
</source>

Example: Minimum Alternate Allele Count = 4
Sample1 Sample2 Sample3 Keep/Discard
0|0 1|1 2|2 Keep
0|0 0|1 2|2 Discard, only 3 Alternates (1 Allele1 & 2 Allele 2)
0|0 1|1 1|2 Keep
0|2 1|1 2|2 Keep
2|1 0|1 2|0 Keep

Example: Minimum Alternate Allele Count = 3 & Exclude Sample2 (without the exclusion, all would be kept)
Sample1 Sample2 Sample3 Keep/Discard
0|0 1|1 2|2 Discard, only 2 Alternates (0 Allele1 & 2 Allele 2)
0|0 0|1 2|2 Discard, only 2 Alternates (0 Allele1 & 2 Allele 2)
0|0 1|1 1|2 Discard, only 2 Alternates (1 Allele1 & 1 Allele 2)
0|2 1|1 2|2 Keep
2|1 0|1 2|0 Keep

===== Minimum Minor Allele Count =====
To Discard any records without a minimum number of minor alleles, use:
<source lang="cpp">
VcfFileReader::addDiscardMinMinorAlleleCount(int32_t minMinorAlleleCount, VcfSubsetSamples* subset)
</source>

The <code>minMinorAlleleCount</code> parameter is the minimum number of minor alleles found in the specified subset (if specified) in order for the record to be kept.

The <code>VcfSubsetSamples* subset</code> parameter is a pointer to the subset of samples that you want to include when counting the number of alleles. If all samples that are read/kept are to be included, NULL should be passed in.

See [[#Handling a Subset of Samples|Handling a Subset of Samples]] for how to use <code>VcfSubsetSamples</code>.

Use the following method to remove the DiscardMinMinorAlleleCount rule:
<source lang="cpp">
VcfFileReader::rmDiscardMinMinorAlleleCount()
</source>

Example: Minimum Minor Allele Count = 2
Sample1 Sample2 Sample3 Keep/Discard
0|0 1|1 2|2 Keep
0|0 0|1 2|2 Discard, only 1 Allele1
0|0 1|1 1|2 Discard, only 1 Allele2
0|2 1|1 2|2 Discard, only 1 Allele0
2|1 0|1 2|0 Keep

Example: Minimum Minor Allele Count = 1 & Exclude Sample2 (without the exclusion, all would be kept)
Sample1 Sample2 Sample3 Keep/Discard
0|0 1|1 2|2 Discard, 0 Allele1
0|0 0|1 2|2 Discard, 0 Allele1
0|0 1|1 1|2 Keep
0|2 1|1 2|2 Discard, 0 Allele1
2|1 0|1 2|0 Keep

==== Read only Certain Sections of the File / Using a VCF Index (TABIX) File ====
The VCF library has the capability of using a tabix file to allow random access in the VCF file.

To use the INDEX File to read sections of the bam file:
# Open the VCF file using: <code>VcfFileReader::open</code>
# Open the VCF Index file using: <code>VcfFileReader::readVcfIndex</code>
# Set the Section to be read: <code>VcfFileReader::setReadSection</code> or <code>VcfFileReader::set1BasedReadSection</code>
# Read records from the section: <code>VcfFileReader::readRecord</code>

<code>VcfFileReader::readVcfIndex</code> returns true if the index file was successfully read and false if not.
<source lang="cpp">
bool readVcfIndex(const char * filename);
bool readVcfIndex();
</source>
If you specify a filename, it looks for the index file with that path/name. If you do not specify the index file name, it will attempt to find it. First it searches for a file with your VCF FileName + a ".tbi" extension. If that isn't found, it removes the ".vcf" if found in your VCF FileName and looks for a file with that name and a ".tbi" extension.

After opening both the VCF File and the Index file, set the section that should be read next using <code>setReadSection</code> or <code>set1BasedReadSection</code>.
<source lang="cpp">
bool setReadSection(const char* chromName);
bool set1BasedReadSection(const char* chromName,
int32_t start, int32_t end,
bool overlap = false);
</source>
Both methods currently always return true, since in the current implementation this can't fail. It currently allows the VcfIndex file to be read after this call.

The read section is reset when a new file is opened, so do not set this prior to opening the file.

<code>setReadSection</code> will set the code to read the entire specified chromosome. On the first <code>readRecord</code> call it will jump to the beginning of the chromosome. It will continue to read the chromosome for each <code>readRecord</code> call made until it has read the entire chromosome. Once the whole chromosome has been read, <code>readRecord</code> will return false until a new read section is specified.

<code>set1BasedReadSection</code> will set the code to read the specified chromosome starting at the specified 1-based position up to, but not including, the specified 1-based end position. On the first <code>readRecord</code> call it will jump the file to the beginning of this section. It will continue to read the section for each <code>readRecord</code> call made until it has read the entire section. Once the entire section has been read, <code>readRecord</code> will return false until a new read section is specified.

<code>set1BasedReadSection</code> takes an optional parameter, <code>overlap</code>. It is defaulted to false which means that only records that start within the specified region will be returned by <code>readRecord</code>. If <code>overlap</code> is set to true, <code>readRecord</code> will also return records that start prior to the specified region, but whose deletions occur in the specified region.

Implementation NOTE: internally it may read extra records prior to the section, but <code>readRecord</code> will keep reading and will not return until it finds a record in the section or until the section has been passed (no records in the section).

<source lang="cpp">
// Open the vcf file for reading.
VcfFileReader reader;
VcfHeader header;
VcfRecord record
reader.open("vcfFileName.vcf", header);

// Index File name is "vcfFileName.vcf.tbi" or "vcfFileName.tbi"
reader.readVcfIndex();

reader.set1BasedReadSection("1", 16384, 32767);
while(reader.readRecord(record))
{
// Process a record from this section.
}

// Set another section, chromosome "X"
reader.setReadSection("X");
while(reader.readRecord(record))
{
// Process record on chromosome X
}
</source>

Alternatively, the index file name can be specified.
<source lang="cpp">
// Index File name is specified
reader.readVcfIndex("myIndex.tbi");
</source>

==== Check if EOF was reached ====
To see if the End Of the File has been reached, use <code>isEOF()</code>.
<source lang="cpp">
reader.isEOF();
</source>
True is returned if the End of the File has been reached, false if not.

==== See How Many Records Were Read ====
To see the total number of records in the file that were read, including any that were discarded, use <code>getNumRecords</code>.

To see only the records in the file that were kept (not discarded) use <code>getNumKeptRecords</code>.
<source lang="cpp">
std::cout << reader.getNumRecords() << " were found in the file\n."
<< "but only " << reader.getNumKeptRecords() << " were processed\n.";
</source>

==== Closing the VCF File ====
<code>close</code> closes the opened VCF file.

<source lang="cpp">
reader.close();
</source>

=== <code>VcfRecord</code> ===
VCF records are read from VCF files using <code>VcfFileReader::readRecord</code>.

Once you have the record, use methods from the <code>VcfRecord</code> to extract the desired information.

==== Extracting Data from the Record ====
Do not delete any of the returned pointers. They are references to strings within the record and will not be valid when the record changes or is deleted.

* Position Information
** <code>getChromStr()</code> returns const char* to this record's chromosome string
** <code>get1BasedPosition()</code> returns an integer containing the 1-based position of the record
* Identifiers
** <code>getIDStr()</code> returns const char* to this record's ID. It includes the entire list of unique identifiers delimited by ';'
* Information on Bases
** <code>getRefStr()</code> returns const char* to the entire reference string
** <code>getAltStr()</code> returns const char* to the entire alternate string (alternates seperated by ',')
** <code>getNumAlts()</code> returns an integer containing the number of alternates in the alt string.
* Quality Information
** <code>getQual()</code> returns a float containing the phred quality score
** <code>getQualStr()</code> returns a string containing the phred quality score
* Filter Information - see [[#VcfRecordFilter|VcfRecordFilter]] for more information on extracting information
** <code>getFilter()</code> returns a reference to the Info field information, a VcfRecordInfo object reference
** <code>passedAllFilters()</code> returns true if the FILTER field indicates <code>PASS</code>, false if not.
* Aditional Information
** <code>getInfo()</code> returns a reference to the filter information, a VcfRecordFilter object reference
* Genotype Information
** <code>getGenotypeInfo()</code> returns a reference to the information from the genotype fields, a VcfRecordGenotype object reference
** <code>allPhased()</code> returns true if all the samples are phased and none unphased and false if any are not phased
** <code>allUnphased()</code> returns true if all the samples are unphased and none phased and false if any are not unphased
** <code>hasAllGenotypeAlleles()</code> returns true if all the samples have all the genotype alleles specified and false if any are missing or the GT field is missing
** <code>getNumSamples()</code> returns the number of samples (that are kept)
** <code>getNumGTs(int index)</code> returns the number of GTs for the specified sample index (starts at 0)
** <code>getGT(int sampleNum, unsigned int gtIndex)</code> returns the integer GT value for the specified sampleNum and GT index (both start at 0). A GT of VcfGenotypeSample::INVALID_GT is returned if the sampleNum or GT index is out of range. A GT of VcfGenotypeSample::MISSING_GT if the GT is '.' (This method is also found in VcfRecordGenotype.

=== VcfRecordFilter ===
<code>VcfRecords</code> contain the data from the <code>INFO</code> field in a <code>VcfRecordFilter</code> object.

== Handling a Subset of Samples ==

When reading a file if you only want to process/keep a subset of samples, use [[#Reading a Subset of Samples|Reading a Subset of Samples]]. When that method is used, only the specified samples are stored. Any further processing will only be on those samples.

Some methods allow the user to specify a subset of samples to operate on. The subset specified when reading the VCF file, if any, is automatically applied since only those samples were stored. If a different/additional subset needs to be applied for other processing, you can use the <code>VcfSubsetSamples</code> class.

To setup a VcfSubsetSamples object, pass the already set VCF header to:
<source lang="cpp">
void VcfSubsetSamples::init(const VcfHeader& header, bool include)
</source>
Set the <code>include</code> parameter to:
* true if all samples should be included except any that are specified as excluded.
* false if all samples should be excluded except any that are specified as included.

NOTE: the header is not modified to add/remove any samples.

To mark a specific sample as excluded use:
<source lang="cpp">
bool VcfSubsetSamples::addExcludeSample(const char* sampleName);
</source>
To mark a specific sample as included use:
<source lang="cpp">
bool VcfSubsetSamples::addIncludeSample(const char* sampleName);
</source>

File:VcfReader.v1.tar

2011-02-04T15:50:34Z

Youna: uploaded a new version of "File:VcfReader.v1.tar": 1. fixed the bug of IDfile. 2. included the functionality to subset variants of a given positions.

2011-02-01T22:10:53Z

Youna:

RvTests

2011-02-01T19:29:43Z

Youna:

[[Category:Software]]
= Overview =
A few rare variants tests (Li-Leal's CMC and Madsen-Browning's weighted method) are implemented in the logisitc regression framework using C++. Please contact Youna Hu (youna@umich.edu) for comments, suggestions or questions.

The source code is located at [[File:RV3Tests.v1.tar]]

You can just download the tar file, extract it and go to the RV3Test.v1 folder to type make all to compile the code, the binary file will then be in the exectuables folder.

= Example =

See a detailed [[example]] here.

= Syntax =

This software uses command line interface as follows

RARE VARIANT ANALYSIS OPTIONS:
GENOTYPE : --genofile [pos.012],
--geneList [outGeneSorted.txt], --cutoff [0.010],
--collapseChoice [or]
PHENOTYPE : --phenofile [LDL.y.ID]
COVARIATES : --covConsider, --covfile [covFile.ID.2.txt]
PERMUTATION : --nPermute [10], --PermutationSeed [1]
GENE LEVEL TEST RESULT : --geneGlobalTestOut [globalPermuteSummary.txt],
--geneTestpvalueFile [geneTestPvalues.txt]

;GENOTYPE

;--genofile: A genotype 012 matrix (.012 is the file) This file can be prepared by using the prepare012s
source code [[File:vcfReader.v1.tar]]
Again, extract the tar file and then go into the directory to type make all to compile the code,
you then will find binary file in the executables folder.

Note: You should use Yanming's vcf annotation [http://genome.sph.umich.edu/wiki/VcfCodingSnps] on your vcf file first
to output a annotated vcf file. You SHOULD keep the log file from the annotation, which will be used to create the gene list.

Data File PREPARATION
Input files : --vcf [LDL.test.vcf], --log [], --IDfile []
Subsetting choices : --All
Output files : --outputPrefix [subsetGeno],
--outputGeneList [LDL.geneList.txt]
--vcf: Input vcf file
--log: This is the log file from Yanming's annotation output, we use this log to obtain the gene list
--IDfile specifies a file with one column of subject IDs to subject from the vcf file.
If it is not specified, then all subjects are included for the format conversion.
--All: specifies 1 to include all variants and 0 to include only nonsyn and stop annotated variants.
-- outputPrefix: Specify the prefix for the four output files which will be used in rvTests
*.012: A genotype matrix with subjects as rows and variant sites as columns.
*.012.pos: Chromosome and position numbers.
*.012.indv: Subject IDs.
*.012.frq: The frequency of the included variants.
--outputGeneList: Specify a file to store the gene list which will be used in rvTest.
The list file looks like this
1 OR4F5 69090 70008
1 SAMD11 860529 871276
1 NOC2L 879583 893918
1 KLHL17 895966 901095
1 PLEKHN1 901876 910482
1 C1orf170 910578 912021

;--geneList: This file is an output from prepare012s using the option --outputGeneList with columns as chromosome number, gene Name, start position, end position. There should be no header for this file.

THE CHROMOSOME NUMBERS SHOULD BE NUMERICS!!!! 1 - chromosome 1, DO NOT USE chr1.

;--cutoff: This is the minor allele frequency, you can specify it as 0.01, 0.05 or etc.
;--collapseChoice: Specify one of {or,sum,wt}. or: Li-Leal's CMC test, sum: Use the number of rare variants for each subject as the score, wt: Madeson-Browning's weighted rare variant score.

;PHENOTYPE
;--phenofile: A file where the first column is subject ID and the second column is phenotype (0 or 1).

;COVARIATES
;--covConsider: Default = 0, no covariate is considered. 1. covariate is considered.
;--covfile: Covariate file with the first column as subject ID and the other columns are covariates needed to be considered in the model.

;PERMUTATION
;--nPermute: Number of permutation for the evaluation of p values.
;-- PermutationSeed: Default = 1. Can be changed to other numbers too.

;GENE LEVEL TEST RESULT:
;--geneGlobalTestOut: This file stores the 5% and 95% quantiles of the p values for all the genes at each permutation
;--geneTestPvalueFile: This file gives you the gene name, number of rare variants, count of variants in case/control and p values from the RV test specified by collapseChoice.

RvTests

2011-02-01T19:29:29Z

Youna:

[[Category:Software]]
= Overview =
A few rare variants tests (Li-Leal's CMC and Madsen-Browning's weighted method) are implemented in the logisitc regression framework using C++. Please contact Youna Hu (youna@umich.edu) for comments, suggestion or questions.

The source code is located at [[File:RV3Tests.v1.tar]]

You can just download the tar file, extract it and go to the RV3Test.v1 folder to type make all to compile the code, the binary file will then be in the exectuables folder.

= Example =

See a detailed [[example]] here.

= Syntax =

This software uses command line interface as follows

RARE VARIANT ANALYSIS OPTIONS:
GENOTYPE : --genofile [pos.012],
--geneList [outGeneSorted.txt], --cutoff [0.010],
--collapseChoice [or]
PHENOTYPE : --phenofile [LDL.y.ID]
COVARIATES : --covConsider, --covfile [covFile.ID.2.txt]
PERMUTATION : --nPermute [10], --PermutationSeed [1]
GENE LEVEL TEST RESULT : --geneGlobalTestOut [globalPermuteSummary.txt],
--geneTestpvalueFile [geneTestPvalues.txt]

;GENOTYPE

;--genofile: A genotype 012 matrix (.012 is the file) This file can be prepared by using the prepare012s
source code [[File:vcfReader.v1.tar]]
Again, extract the tar file and then go into the directory to type make all to compile the code,
you then will find binary file in the executables folder.

Note: You should use Yanming's vcf annotation [http://genome.sph.umich.edu/wiki/VcfCodingSnps] on your vcf file first
to output a annotated vcf file. You SHOULD keep the log file from the annotation, which will be used to create the gene list.

Data File PREPARATION
Input files : --vcf [LDL.test.vcf], --log [], --IDfile []
Subsetting choices : --All
Output files : --outputPrefix [subsetGeno],
--outputGeneList [LDL.geneList.txt]
--vcf: Input vcf file
--log: This is the log file from Yanming's annotation output, we use this log to obtain the gene list
--IDfile specifies a file with one column of subject IDs to subject from the vcf file.
If it is not specified, then all subjects are included for the format conversion.
--All: specifies 1 to include all variants and 0 to include only nonsyn and stop annotated variants.
-- outputPrefix: Specify the prefix for the four output files which will be used in rvTests
*.012: A genotype matrix with subjects as rows and variant sites as columns.
*.012.pos: Chromosome and position numbers.
*.012.indv: Subject IDs.
*.012.frq: The frequency of the included variants.
--outputGeneList: Specify a file to store the gene list which will be used in rvTest.
The list file looks like this
1 OR4F5 69090 70008
1 SAMD11 860529 871276
1 NOC2L 879583 893918
1 KLHL17 895966 901095
1 PLEKHN1 901876 910482
1 C1orf170 910578 912021

;--geneList: This file is an output from prepare012s using the option --outputGeneList with columns as chromosome number, gene Name, start position, end position. There should be no header for this file.

THE CHROMOSOME NUMBERS SHOULD BE NUMERICS!!!! 1 - chromosome 1, DO NOT USE chr1.

;--cutoff: This is the minor allele frequency, you can specify it as 0.01, 0.05 or etc.
;--collapseChoice: Specify one of {or,sum,wt}. or: Li-Leal's CMC test, sum: Use the number of rare variants for each subject as the score, wt: Madeson-Browning's weighted rare variant score.

;PHENOTYPE
;--phenofile: A file where the first column is subject ID and the second column is phenotype (0 or 1).

;COVARIATES
;--covConsider: Default = 0, no covariate is considered. 1. covariate is considered.
;--covfile: Covariate file with the first column as subject ID and the other columns are covariates needed to be considered in the model.

;PERMUTATION
;--nPermute: Number of permutation for the evaluation of p values.
;-- PermutationSeed: Default = 1. Can be changed to other numbers too.

;GENE LEVEL TEST RESULT:
;--geneGlobalTestOut: This file stores the 5% and 95% quantiles of the p values for all the genes at each permutation
;--geneTestPvalueFile: This file gives you the gene name, number of rare variants, count of variants in case/control and p values from the RV test specified by collapseChoice.

File:VcfReader.v1.tar

2011-01-31T17:18:02Z

Youna: uploaded a new version of "File:VcfReader.v1.tar"

2011-01-30T20:00:42Z

Youna:

[[Category:Software]]
= Overview =
A few rare variants tests (Li-Leal's CMC and Madsen-Browning's weighted method) are implemented in the logisitc regression framework using C++.

The source code is located at [[File:RV3Tests.v1.tar]]

You can just download the tar file, extract it and go to the RV3Test.v1 folder to type make all to compile the code, the binary file will then be in the exectuables folder.

= Example =

See a detailed [[example]] here.

= Syntax =

This software uses command line interface as follows

RARE VARIANT ANALYSIS OPTIONS:
GENOTYPE : --genofile [pos.012],
--geneList [outGeneSorted.txt], --cutoff [0.010],
--collapseChoice [or]
PHENOTYPE : --phenofile [LDL.y.ID]
COVARIATES : --covConsider, --covfile [covFile.ID.2.txt]
PERMUTATION : --nPermute [10], --PermutationSeed [1]
GENE LEVEL TEST RESULT : --geneGlobalTestOut [globalPermuteSummary.txt],
--geneTestpvalueFile [geneTestPvalues.txt]

;GENOTYPE

;--genofile: A genotype 012 matrix (.012 is the file) This file can be prepared by using the prepare012s
source code [[File:vcfReader.v1.tar]]
Again, extract the tar file and then go into the directory to type make all to compile the code, you then will find binary file in the executables folder.

Note: If you going to analyze nonsynonymous and stop annotated variants,
you should use Yanming's vcf annotation [http://genome.sph.umich.edu/wiki/VcfCodingSnps] on the vcf file.

Data File PREPARATION
Input files : --vcf [LDL.test.vcf], --log [], --IDfile []
Subsetting choices : --All
Output files : --outputPrefix [subsetGeno],
--outputGeneList [LDL.geneList.txt]
--vcf: Input vcf file
--log: This is the log file from Yanming's annotation output, we use this log to obtain the gene list
--IDfile specifies a file with one column of subject IDs to subject from the vcf file.
If it is not specified, then all subjects are included for the format conversion.
--All: specifies 1 to include all variants and 0 to include only nonsyn and stop annotated variants.
-- outputPrefix: Specify the prefix for the four output files which will be used in rvTests
*.012: A genotype matrix with subjects as rows and variant sites as columns.
*.012.pos: Chromosome and position numbers.
*.012.indv: Subject IDs.
*.012.frq: The frequency of the included variants.
--outputGeneList: Specify a file to store the gene list which will be used in rvTest.
The list file looks like this
1 OR4F5 69090 70008
1 SAMD11 860529 871276
1 NOC2L 879583 893918
1 KLHL17 895966 901095
1 PLEKHN1 901876 910482
1 C1orf170 910578 912021

;--geneList: This file is an output from prepare012s using the option --outputGeneList with columns as chromosome number, gene Name, start position, end position. There should be no header for this file.

THE CHROMOSOME NUMBERS SHOULD BE NUMERICS!!!! 1 - chromosome 1, DO NOT USE chr1.

;--cutoff: This is the minor allele frequency, you can specify it as 0.01, 0.05 or etc.
;--collapseChoice: Specify one of {or,sum,wt}. or: Li-Leal's CMC test, sum: Use the number of rare variants for each subject as the score, wt: Madeson-Browning's weighted rare variant score.

;PHENOTYPE
;--phenofile: A file where the first column is subject ID and the second column is phenotype (0 or 1).

;COVARIATES
;--covConsider: Default = 0, no covariate is considered. 1. covariate is considered.
;--covfile: Covariate file with the first column as subject ID and the other columns are covariates needed to be considered in the model.

;PERMUTATION
;--nPermute: Number of permutation for the evaluation of p values.
;-- PermutationSeed: Default = 1. Can be changed to other numbers too.

;GENE LEVEL TEST RESULT:
;--geneGlobalTestOut: This file stores the 5% and 95% quantiles of the p values for all the genes at each permutation
;--geneTestPvalueFile: This file gives you the gene name, number of rare variants, count of variants in case/control and p values from the RV test specified by collapseChoice.