Dinh/Dinh 2012/NOTES/2012-3-22
Jump to navigation
Jump to search
BSPP data trimming and cleaning rules using 48 African (UPenn) as an example
Methylation
- Script to generate methylation matrix: File:AllBED2Matrix DD.txt
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
- Once a cut-off is determined, use awk and cut to get the correct matrix file
awk '{if($46 > 0.2363) print $0}' autosomes.clean.UPenn45.hg19.min40.methylMatrix | cut -f1-46 > pca.UPenn44.hg19.min40.methylMatrix
BSPP SNP
- Script to generate SNP matrix (.tped/.tfam): File:GetTPED.txt
getTPED.pl UPenn44.filteredSNP134 < list_snp44
- First argument is the file name suffix
- list_snp44 : simply list path to all included SNP files.
- NOTE:
plink doesn't link it when a locus have a reported number of alleles which is greater than 2. This will result in an error message during conversion. We don't have this problem for the array SNP, probably because they have been pre-filtered but this is a problem for BSPP SNPs. To filter for these errors, in the script getTPED.pl, 1) Count all alleles 2) If number of alleles is 4, then don't report this locus in the matrix 3) If number of alleles is 3, find the lowest frequency allele (if there is a tie, pick the last one), and mask any SNP calls containing this allele (convert to NAs or "0 0")
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")