Dinh/Dinh 2012/NOTES/2012-3-22: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
mNo edit summary
>Dinh
mNo edit summary
 
(3 intermediate revisions by the same user not shown)
Line 1: Line 1:
=BSPP data trimming and cleaning rules using 48 African (UPenn) as an example=
=BSPP data trimming and cleaning procedure using 48 African (UPenn) as an example=


==Methylation==
==Methylation==
Line 37: Line 37:
  getBestSites.pl 49 46 < clean.UPenn44.hg19.min30.methylMatrix > tmp
  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  
* 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==
==BSPP SNP==
* Script to generate SNP matrix (.tped/.tfam): [[File:getTPED.txt]]
* Script to generate SNP matrix (.tped/.tfam): [[File:getTPED.txt]]
Line 43: Line 44:
* First argument is the file name suffix
* First argument is the file name suffix
* '''list_snp44''' : simply list path to all included SNP files.
* '''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===
===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):  
* Filter matrix to allow maximum 25% missing values (25% for QTL, may require less for clustering):  
Line 54: Line 60:
   (plink.mdist will be generated, rename the file:)
   (plink.mdist will be generated, rename the file:)
   mv plink.mdist UPenn44_IBD
   mv plink.mdist UPenn44_IBD
* Note that this matrix does not have headers for both rows and columns, so when clustering using R
* Note that this matrix does not have headers for both rows and columns, so when clustering using R, we need to add-in this information
   A = read.table("*.tfam",F)
   A = read.table("*.tfam",F)
   B = read.table("UPenn44_IBD",F)
   B = read.table("UPenn44_IBD",F)

Latest revision as of 00:29, 29 March 2012

BSPP data trimming and cleaning procedure using 48 African (UPenn) as an example[edit]

Methylation[edit]

 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[edit]

  • 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[edit]

  • 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[edit]

  • 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[edit]

  • 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[edit]

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[edit]

  • 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, we need to add-in this information
 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")