Since R can automatically recognize and read compressed files with standard name extensions like .gz or .bzip2, this holds as well for the input files that can be used with this package.
Loading data files: the function data2haplohh
The data2haplohh() function converts data files into an R object of class haplohh, forming the basis for subsequent analysis. The following options are available to recode alleles and/or to filter markers and haplotypes:
1. Allele coding options:
The parameter allele_coding has four options: three of them expect the user to provide a coding, while the fourth option overwrites any user defined coding.
The package accepts two different integer codings within the haplotype file which can be specified as "12" or "01", respectively:
| Missing value |
0 |
NA or . |
| Ancestral allele |
1 |
0 |
| Derived allele 1 |
2 |
1 |
| Derived allele 2 |
3 |
2 |
| … |
… |
… |
The latter one is also the coding used internally by our package.
If the alleles are given in an un-coded form (e.g. as nucleotides), the alleles can be coded for each marker separately with help of the marker information file. This option is activated by allele_coding = "map". If an allele of a marker is found in the fourth column of the marker information file, it is re-coded as 0 (ancestral) and if it is found in the fifth column it is re-coded as 1 (derived). The fifth column may contain several alleles, separated by commas, like in the variant call format. In this case, the internal coding extends to higher integers. If an allele present in the haplotype file is not found neither in the fourth nor the fifth column of the marker information file, it is counted as NA (missing value).
If the alleles are not polarized, i.e. it is unknown which allele is ancestral and which derived, the option "none" should be selected. The alleles will then be coded in alpha-numeric order for each marker; the first becomes 0, the second 1, etc. Missing values must be given as NA or ‘.’ while any other number, character or string will be counted as an allele. Options of subsequent analysis functions have to be set accordingly (see section 7.6).
2. Filtering options
2.1 Discard haplotypes with a high amount of missing data. If haplotypes contain missing data (presumably a rare case since most phasing programs allow imputing missing genotypes), it is possible to discard those with a fraction of less than min_perc_geno.hap of markers genotyped. Its default is NA, hence no restriction.
2.2. Discard markers with a high amount of missing data. It is possible to discard individual markers genotyped on a fraction of less than min_perc_geno.mrk of haplotypes. By default min_perc_geno.mrk = 100 meaning that only fully genotyped markers are retained. A value of NA or zero is not allowed since rehh cannot handle markers with no data.
2.3 Discard markers with a low Minor Allele Frequency. It is possible to discard markers with a Minor Allele Frequency (MAF) equal to or below min_maf. Its default is NA meaning no constraint. Setting this value to zero eliminates monomorphic markers. Note that since low-frequency alleles usually abound, any positive value may exclude a substantial fraction of the data.
The arguments min_perc_geno.hap, min_perc_geno.mrk and min_maf are evaluated in this order.
More details about input options can be found in the documentation using the command:
Example 1: reading haplotype file in standard format
The following command converts the haplotype input file bta12_cgu.hap (in “standard” haplotype format) and marker information input file map.inp to an haplohh object named hh. Because the map file contains information for markers mapping to multiple chromosomes we have to specify by chr.name = 12 that the haplotype input file is about chromosome 12. Alleles are given as nucleotides in the haplotype file and coded with help of the map file by setting allele_coding = "map".
hh <- data2haplohh(hap_file = "bta12_cgu.hap",
map_file = "map.inp",
chr.name = 12,
allele_coding = "map")
> * Reading input file(s) *
> Map info: 1424 markers declared for chromosome 12 .
> Haplotype input file in standard format assumed.
> Alleles are being recoded according to fourth and fifth column of map file.
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Data consists of 280 haplotypes and 1424 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 27 1397
If no value is specified for chr.name and more than one chromosome is detected in the map file, an error is produced:
hh <- data2haplohh(hap_file = "bta12_cgu.hap",
map_file = "map.inp",
allele_coding = "map")
> * Reading input file(s) *
> More than one chromosome name in file: map.inp .
> Here is a list of chromosome names found:
> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
> Error: Please specify a chromosome name. Conversion stopped.
Furthermore, the following message is prompted if the number of markers for the chromosome in the info file does not correspond to the one in the haplotype file (for instance when a wrong chromosome is specified):
hh <- data2haplohh(hap_file = "bta12_cgu.hap",
map_file = "map.inp",
chr.name = 18,
allele_coding = "map")
> * Reading input file(s) *
> Map info: 1123 markers declared for chromosome 18 .
> Haplotype input file in standard format assumed.
> Error: The number of markers in the haplotype file (1424) is not equal
> to the number of markers in the map file (1123).
> Conversion stopped.
Example 2: reading haplotype file in transposed format (SHAPEIT2–like)
If the haplotype input file is in “transposed” format (like bta12_cgu.thap), the option haplotype.in.columns has to be set to TRUE while all other parameters remain unaffected with respect to example 1. This is the only format which is not recognized automatically, but has to be explicitly declared by the user.
hh <- data2haplohh(hap_file = "bta12_cgu.thap",
map_file = "map.inp",
chr.name = 12,
allele_coding = "map",
haplotype.in.columns = TRUE)
> * Reading input file(s) *
> Map info: 1424 markers declared for chromosome 12 .
> Haplotype input file in transposed format assumed.
> Alleles are being recoded according to fourth and fifth column of map file.
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Data consists of 280 haplotypes and 1424 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 27 1397
Example 3: reading haplotype file in fastPHASE output format
In this example, we use the fastPHASE output bta12_hapguess_switch.out as haplotype file. This format is recognized automatically. As before, we need to set chr.name = 12.
Because haplotypes originated here from several populations (the fastPHASE option -u was used), it is necessary to select the population of interest (in our case “CGU”)
by its corresponding number in the file, hence popsel = 7.
hh <- data2haplohh(hap_file = "bta12_hapguess_switch.out",
map_file = "map.inp",
chr.name = 12,
popsel = 7,
allele_coding = "map")
> * Reading input file(s) *
> Map info: 1424 markers declared for chromosome 12 .
> Haplotype input file in fastPHASE format assumed.
> Haplotypes in the fastPHASE output file originate from 8 populations.
> Alleles are being recoded according to fourth and fifth column of map file.
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Data consists of 280 haplotypes and 1424 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 27 1397
If no value is specified for the popsel argument and more than one population is detected in the fastPHASE output file, an error is produced and the available population numbers printed:
hh <- data2haplohh(hap_file = "bta12_hapguess_switch.out",
map_file = "map.inp",
chr.name = 12,
allele_coding = "map")
> * Reading input file(s) *
> Map info: 1424 markers declared for chromosome 12 .
> Haplotype input file in fastPHASE format assumed.
> Haplotypes in the fastPHASE output file originate from 8 populations.
> Error: Please specify by 'popsel' one of the following population numbers:
> 1 2 3 4 5 6 7 8
> Conversion stopped.
Example 4: reading vcf files
The function data2haplohh() checks automatically whether the specified haplotype input file is in variant call format (vcf). If this is the case, the parameters map_file and allele_coding are ignored. By default, the function tries to polarize (see section 7.6) the alleles of each marker using the ancestral allele, expected to be given by key AA of the INFO field.
Ancestral alleles are sometimes marked by upper case as “high confident” and by lower case as “low confident”. The default setting capitalize_AA = TRUE lifts this distinction before polarization.
If the AA key is absent, the option polarize_vcf should be set to FALSE and the allele coding of the vcf file is directly used as internal coding.
If there is data for more than one chromosome in the file, the chromosome of interest has to be specified by chr.name. Since always the whole file is read in, it is advisable to split large data sets into chromosome-specific files.
In order to process vcf files, the package vcfR or the package
data.table (which in turn needs R.utils to read compressed files) have to be installed. The parameter vcf_reader has to be set to either "vcfR" or "data.table".
In the file bta12_cgu.vcf.gz the ancestral allele was set as reference and hence no further polarizing is necessary.
hh <- data2haplohh(hap_file = "bta12_cgu.vcf.gz",
polarize_vcf = FALSE,
vcf_reader = "data.table")
> * Reading input file(s) *
> Using package 'data.table' to read vcf.
> Extracting map information.
> Extracting haplotypes.
> Number of individuals which are
> Haploid Diploid Triploid, ... :
> 1 2
> 0 140
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Data consists of 280 haplotypes and 1424 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 27 1397
Example 5: reading ms output
The function data2haplohh() automatically checks whether the haplotype file is in the output format of the simulation program ms (Hudson 2002). If this is the case, the parameters map_file and allele_coding are ignored. If the file contains several ‘runs’ (as referred to by the parameter nrep of ms), it is necessary to specify the number of the run in option chr.name. Note that always the whole file is read, so that it might be advisable to spread large simulations over separate files.
One argument of the data2haplohh function is specifically dedicated to ms output, although it works with other formats as well: ms gives chromosomal positions as fractions of the interval [0,1] and in order to obtain more realistic values, these positions can be multiplied by a factor, set by position_scaling_factor.
Note that ms output can contain multiple markers with the same (rounded) position, which rehh does not accept. In this case the numerical precision for chromosomal positions in the ms output should be increased (option -p of ms, option -oformat of msms).
Setting remove_multiple_markers to TRUE entails that from consecutive markers with the same position only the first one is retained and a warning containing the number of removed markers is printed. Note that this effectively transforms the “infinite sites model” used for simulations by ms into a “finite sites model”.
In order to read the ms format, the package gap has to be installed.
hh <- data2haplohh(hap_file = "ms.out",
chr.name = 2,
position_scaling_factor = 1000)
> * Reading input file(s) *
> Input file in 'ms' output format assumed.
>
> Extracting positions.
> Extracting haplotypes.
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Data consists of 10 haplotypes and 39 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 0 39
Subset data
Sometimes it might be of interest to restrict the computation of statistics to subsets of haplotypes or markers. This can be done with help of the function subset() and arguments select.hap and select.mrk. After subsetting, the same filters are applied as during data input from files in order to ensure that no marker has exclusively missing values and, optionally, to exclude markers which are monomorphic within the subset by setting min_maf to zero.
For example, to restrict the output of ms (see section 2.3.5) to the first five haplotypes and filter out monomorphic sites, we apply
hh_subset = subset(hh, select.hap = 1:5, min_maf = 0)
> * Subset haplotypes *
> 5 haplotypes discarded.
> 5 haplotypes remaining.
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Discard markers with Minor Allele Frequency equal to or below 0 .
> 30 markers discarded.
> 9 markers remaining.
> Data consists of 5 haplotypes and 9 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 0 9
while the following command omits the first marker:
hh_subset = subset(hh, select.mrk = -1)
> * Subset markers *
> 1 markers discarded.
> 38 markers remaining.
> * Filtering data *
> Discard markers genotyped on less than 100 % of haplotypes.
> No marker discarded.
> Data consists of 10 haplotypes and 38 markers.
> Number of mono-, bi-, multi-allelic markers:
> 1 2
> 0 38
