Dinh 2011/NOTES/2011-11-18
Jump to navigation
Jump to search
[edit]
- I did the mQTL on the 30 samples which we have genotyping array data for (geno_GA.txt)
- Also SOM analysis on these 30 samples
- All 30 are diseased patients ... Not useful for finding bio-markers
Generating tped and tfam and methylation matrix files[edit]
- TFAM
awk '{print $2}' geno_GA.txt | sort -u > geno_GA.tfam (removed unwanted lines first) awk '{print sum" "$2" 0 0 0 -9"; sum++}' geno_GA.tfam > tmp vim tmp # add in number 30 for the first row. sort -k1,1n tmp > geno_GA.tfam
- getChrPos.pl for rs sites
./getChrPos.pl geno_GA.txt < ../HsRef/Hg19/snp134_snv.txt > annotated.geno_GA.txt &
- getTPED_11172011.pl
./getTPED_11172011.pl geno_GA.tfam < annotated.geno_GA.txt > geno_GA.tped &
- Modified allBED2Matrix to get methylation matrix:
The script was modified to read only 30 specified samples nohup ./allBED2Matrix.pl AMD-BED 27 1 0.1 > geno_GA_30.methylMatrix_min27_std0.1 &
- Split CpGs file and run mQTL:
split -l 10000 geno_GA_30* geno_methylMatrix_frag head -1 geno_GA_30* > header for f in geno_methylMatrix_frag* do cat header $f > tmp ./createPhenoFile.pl geno_GA.tfam tmp > $f done; mkdir phenoAB mkdir phenoAC mkdir phenoAD mkdir phenoAE mkdir phenoAF mkdir phenoAG mkdir phenoAH mkdir phenoAI mkdir phenoAJ mkdir phenoAK mkdir phenoAL mkdir phenoAM mv geno_methylMatrix_fragaa phenoAA/matrix.txt mv geno_methylMatrix_fragab phenoAB/matrix.txt mv geno_methylMatrix_fragac phenoAC/matrix.txt mv geno_methylMatrix_fragad phenoAD/matrix.txt mv geno_methylMatrix_fragae phenoAE/matrix.txt mv geno_methylMatrix_fragaf phenoAF/matrix.txt mv geno_methylMatrix_fragag phenoAG/matrix.txt mv geno_methylMatrix_fragah phenoAH/matrix.txt mv geno_methylMatrix_fragai phenoAI/matrix.txt mv geno_methylMatrix_fragaj phenoAJ/matrix.txt mv geno_methylMatrix_fragak phenoAK/matrix.txt mv geno_methylMatrix_fragal phenoAL/matrix.txt mv geno_methylMatrix_fragam phenoAM/matrix.txt for f in phenoA* do cd $f /media/2TB_store2/Dinh_Scratch/scripts/mQTL_2MB_plink.pl ../geno_GA matrix.txt > results.txt & cd .. done;
- I still need to process the result files
- Since this analysis was done on all GA (disease status samples), it might not be worthwhile to continue looking at this.
R package oposSOM analysis[edit]
- Wirth, Loeffler, v.Bergen, Binder: Expression cartography of human tissues using self organizing maps. (BMC Bioinformatics 2011)
library(oposSOM) A = read.table("geno_GA_30.methylMatrix_min27_std0.1", header=TRUE, row.name=1) A$STDEV=NULL A$min_RD=NULL A$mean_RD=NULL B = as.matrix(A) run.pipeline(B, dataset.name="AMD_std0.1_30", dim.som1=20, dim.som2=20)
- Understanding the dim.som1 and dim.som2
- dim.som1 = size of the primary SOM / resolution of the expression profiles: number of tiles equals dim.som1^2.
- dim.som2 = resolution of the second level SOM: number of tiles equals dim.som2^2
- Since the analysis doesn't allow for missing values, can we iteratively switch out samples to find a converging best minimal set?
- To allow for ~4% missing values, I can switch out 2 samples in control and in GA at each iteration.
- The number of ways 2 samples could be selected out of 48: 48!/(48-2)! = 2256. Since, there are control and GA, = 2256^2 = 5,089,536
- Do 1 million random permutation to sample 20% of the solution space? But since we only allow for 4% missing values, there should be very little variations between the solutions. But with probe capture, we only get around 80-90% sites shared between two samples.
- Lets try permuting two solutions each, from 35+35, 37+37, 39+39, 41+41, 43+43, 45+45, 47+47 to see how much overlaps are there between two sets.