Dinh/Dinh 2012/NOTES/2012-3-22

From ZhangLabWiki
Revision as of 21:31, 28 March 2012 by >Dinh
Jump to navigation Jump to search

BSPP data trimming and cleaning rules using 48 African (UPenn) as an example

Methylation

 allBED2Matrix_DD.pl list_bed44 30 1 0 &
  • list_bedXX file have two columns: 1) sampleID, 2) path to BED file
  • 30 is the minimum sample, set strictly to allow maximum 10% missing value per CpG.
  • From our March 16 meeting, we finalized ways to clean data using the simulation results from the 9 technical replicates.

Inclusion/Exclusion rule

  • Take the array of QC (pearson's correlation value) and perform a Grubbs' test
Reject samples with low outlying QC value by p<0.05.
  • Take the array of number of CpGs characterized within each BED file and perform a Grubbs' test
Reject samples with low outlying CpG captured by p<0.05.
  • Rationale: Samples with low QC are problematic because they contain many low-quality sites. We should also reject samples with low CpG captured because this means that the capture did not worked well and there may be many PCR duplicates.

Remove low confidence sites

  • The methylation matrix needs to be cleaned by removing all high technical variability sites:
 Accept site if techRepSTDEV < 0.1 OR p_value > 0.1
  • Since not all CpG sites have been characterized by the 9 technical replicates (only ~200K can be used each time), we need to figure out a way to infer techRevSTDEV and p_value from probe efficiency measures.
  • Script to append p-value from techical replicate: File:Getsite.txt
getsite.pl UPenn44.hg19.min30.methylMatrix < ../Tech_Rep_9/SimulationResults_2172012_header > UPenn44.hg19.min30.methylMatrix.pvalues
  • Command to cleanup:
awk '{if($49 > 0.1 || $50 > 0.1) print $0}' UPenn44.hg19.min30.methylMatrix.pvalues > clean.UPenn44.hg19.min30.methylMatrix

mQTL

  • As a standard requirement, CpGs for mQTL should contain no more than 95% individuals with hypo (<0.2) or hypermethylation (>0.8) since such sites are probably within CpG islands.
  • Script to filter by hypo or hypermethylation: File:FilterbyHyperHypo.txt
filterbyHyperHypo.pl 44 < clean.UPenn44.hg19.min30.methylMatrix | cut -f1-45 > qtl.UPenn44.hg19.min30.methylMatrix
  • First argument is the total number of samples.

PCA

  • As a standard to perform PCA, we need to calculate a subset of sites with low error.
  • Find STDEV = x such that number of sites with STDEV=x in data (A) is greater than the number of sites with STDEV=x in technical replicates (B)
  • When using err ratio (B/A) <= 0.2, this means that 20% of variabilities are contributed by sampling variabilities (maybe some technical variabilities as well).
  • Script to find x: File:GetBestSites.txt
getBestSites.pl 49 46 < clean.UPenn44.hg19.min30.methylMatrix > tmp
  • First argument is the column # (1based) for techRepSTDEV and second is the column# for the data STDEV

BSPP SNP

getTPED.pl UPenn44.filteredSNP134 < list_snp44
  • First argument is the file name suffix
  • list_snp44 : simply list path to all included SNP files.

Use PLINK for trimming SNP matrix and distance calculations

  • Filter matrix to allow maximum 25% missing values (25% for QTL, may require less for clustering):
(suffix name for tped/tfam is UPenn44.filteredSNP134)
 plink-1.07-x86_64/plink --noweb --tfile UPenn44.filteredSNP134 --geno 0.25 --recode --out cleaned
  • Tranpose the clean matrix to tped/tfam again:
 plink-1.07-x86_64/plink --noweb --file cleaned --recode --transpose --out min33.UPenn44.filteredSNP134
  • To get the IBD distance matrix:
 plink-1.07-x86_64/plink --tfile min33.UPenn44.filteredSNP134 --cluster --distance-matrix --noweb
 (plink.mdist will be generated, rename the file:)
 mv plink.mdist UPenn44_IBD
  • Note that this matrix does not have headers for both rows and columns, so when clustering using R
 A = read.table("*.tfam",F)
 B = read.table("UPenn44_IBD",F)
 rownames(B) = A$V2
 colnames(B) = A$V2
 plot(hclust(as.dist(B), method="average"), main="SNP cluster, n=44")