1.1   2008-07-04 check out the intensity histogram of QC probes in CHLA-06-27

Arrays coming out of CHLA on 06-27-2008 have 4 bad ones (mismatch rate >=12%) out of 9 arrays. DNA quantity and call NA rate don't seem to perfectly correlate with mismatch rate. Finally, an ad-hoc badness of intensity histogram, which roughly corresponds to how well the histogram conforms to that of the 107 normalized arrays, predicts the mismatch rate well.

have 3 doubts:

1. the array quality could be bad. magnus said it's expired in April. or it could be something happened in my trip from chicago to here, i put them in a luggage bag, the x-ray machine of airport security could have damaged them ...
2. our calling algorithm needs to adjust to the new intensity distribution of these arrays.
3. scanner alignment is not well.

For 1, we can't really do anything. hope to tell it apart from the image.

For 2, we'll try and see.

For 3, Aaron suggested to look at the boundaries of the raw images.

Eye-balling of the raw images showed the arrays are well-aligned but are not as bright as those chicago arrays. Also bad arrays are also less brighter than the good arrays. This leads to 2 possibilities:

1. DNA quantity is not enough for bad arrays. The numbers in the table were not accurate?
2. DNA hybridization doesn't work well for bad arrays. Washed away?? Check what the intensity histogram of QC probes look like?

Procedures to get QC probe intensity histogram

1. commandline to output QC probe intensity:

   crocea@banyan:~/script/affy/sdk/calvin_files$ for ((i=498;i<507;i++));
   do  echo $i; filename=/Network/Data/250k/db/raw_data/$i\_raw_data.cel;
   ./OutputQCIntensity -i $filename -d \
   ~/script/variation/genotyping/250ksnp/data/atSNPtilx520433_rev2/Full/atSNPtil_geno/LibFiles/atSNPtil_geno.cdf \
   -o /tmp/CHLA-06-27/$i.tsv;
   done
   

2. command line to draw the intensity (internally calling above to get intensity matrix):

   ./src/PlotQCProbeIntensityHistogram.py -i 498 -o /tmp/qc_plots/498.png -m 50000
   

The 9 histograms are in qc-probe-intensity-histogram/ in the order of array-id (Field 'id' in table stock_250k.array_info). There's no difference between bad arrays (id=498, 501, 504, 506) and good arrays. This means:

1. QC probes work fine.
2. DNA hybridization works fine. Otherwise, QC probes should have bad hybridization in bad arrays, which result in lower intensity.
3. The DNA quantity is lower in bad arrays. The sharp peak in the intensity histogram of all probes in bad arrays, which the "badness" measure is based upon, represents noise?

On the other hand, compare the QC probe histogram between CHLA and chicago, although not always, the ones from chicago , peak >20000, are right-shifted (bigger) as to CHLA arrays, peak around 15000, which suggests the overall CHLA array intensity is lower than chicago.

1.2   2009-3-10 CHLA 01-2009 batch

QC for total 429 arrays (ecotype id of 1 array isn't resolved). 104 (=24.2%) has mismatch rate above 10%. a few might be ecotype id mis-linking (might be mislabelled upstream procedures, not necessarily glenda's fault) because the NA (non-calling) rate is pretty low (26 out of 104 has NA rate below 10%). but most of 104 have high NA rates.

391 arrays have jpegs. (430 arrays total though). ~24 have either a curvy line crossing the array or a small smudge, ~13 have large smudge. there're also some arrays substantially less brighter that others.

1.3   2009-4-16 27 mis-labelled arrays among all current arrays

The way to get this is to compare 250k calls of failed arrays against all 149SNP data. True identity is called if the mismatch rate is <=2% and the number of corresponding pairs is >=40. Sometimes, one failed 250k array could be matched to dozens of 149SNP data. Then all of them are listed. so explains why there are 215 of them. They actually correspond to only 27 failed runs.

code:

call_method_id = 3
qc_method_id = 9
output_fname = '/tmp/250k_good_bad_cross_label_arrays.tsv'
reportBadArrays(db_250k, call_method_id, qc_method_id, output_fname)

2.1   Data & Methods

Two sample rank sum test, aka Mann-Whitney (use R's wilcox.test). P-values from the candidate gene list is one sample. P-values from non-candidate genes is another sample. P-value for a gene is defined as the p-value from the most-significant SNP within the gene or nearby (within 50kb). Association between genes and SNPs is defined in table snps_context.

Candidate genes are stored in table candidate_gene_list and candidate_gene_list_type.

2.2   Results

2.3   Conclusions

as far as that flowering time short list is not very bogus (suzi?), a few conclusions are obvious:

1. margarita is not that great.

   One question here is whether Margarita tends to miss the 50 kb window, or simply find more population structure. Either way, I don't view Margarita as one of the main methods. We're basically running it to see if it finds something the others miss.

2. ~170 accession is better than 96 accessions.

3. random forest is unstable. for LD, it's perfect. but not for the other 3 phenotype.

4. kinship matrix generated by 2010 alone is less optimal? K.2010 and K.2010_u2_250K gave same stuff.

3.1   commandlines

submit results from one analysis method but for all phenotypes into db:

for ((i=1;i<=186;i++));
do echo $i;
~/script/variation/src/Results2DB_250k.py -i \
/Network/Data/250k/tmp-bvilhjal/emma_results/Emma_trans_censored_Kinship_250K_i1_2010_$i\_* \
-l 9 -a 6 -e $i -f Emma_trans_censor_ana9_ph$i -u yh -c -p pass**d ;
done;

output rank test p-values for several association results:

./script/variation/src/GeneListRankTest.py -e \
19,20,21,22,571,572,573,574,389,390,391,392,204,205,206,207,943,\
944,1045,1046,1151,1152,1153,1154,1337,1338,1339,1340,1522,1523,1524,1525,1707,1708,1709,1710 \
-l 1 -o ~/Desktop/rank_test_results_pheno1234_noRF.tsv  -u yh

4.1   check missing 192 in 250k

1. thru stock.batch_ecotype to see how many are genotyped in 250k

   select distinct b.ecotypeid from stock.batch_ecotype b, stock.ecotype_duplicate2tg_ecotypeid et, view_qc v where v.ecotype_id=et.tg_ecotypeid and b.batchid=2 and b.ecotypeid=et.ecotypeid and v.call_method_id=3;

2. thru at.magnus_192_vs_accession to see how many are genotyped in 250k

   select distinct a.* from at.magnus_192_vs_accession a, at.accession2tg_ecotypeid ae, call_info c, array_info ar where ar.ecotype_id=ae.ecotype_id and a.accession_id=ae.accession_id and c.method_id=3 and c.array_id=ar.id;

3. show all 192, with not genotyped shown as NULL ecotype_id

   select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar where c.method_id=3 and c.array_id=ar.id) order by ecotype_id, linename;

4. show all 192, the ones not genotyped or genotyped but mismatch_rate>=0.12 having NULL ecotype_id

   select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar, call_qc q where q.call_info_id=c.id and c.method_id=3 and c.array_id=ar.id and q.qc_method_id=9 and q.mismatch_rate<0.12) order by ecotype_id, linename;

5. to get only the ones not genotyped or genotyped but mismatch_rate>=0.12

   select newt.accession_id, ae.nativename, ae.stockparent, ae.ecotype_id from (select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar, call_qc q where q.call_info_id=c.id and c.method_id=3 and c.array_id=ar.id and q.qc_method_id=9 and q.mismatch_rate<0.12) order by ecotype_id, linename) as newt, at.accession2tg_ecotypeid ae where newt.ecotype_id is NULL and newt.accession_id=ae.accession_id;

6. get the Qc info

   select table2.*, q.call_info_id, q.QC_method_id, q.NA_rate, q.no_of_totals, q.mismatch_rate, q.no_of_mismatches, q.no_of_non_NA_pairs from (select newt.accession_id, ae.nativename, ae.stockparent, ae.ecotype_id from (select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar, call_qc q where q.call_info_id=c.id and c.method_id=3 and c.array_id=ar.id and q.qc_method_id=9 and q.mismatch_rate<0.12) order by ecotype_id, linename) as newt, at.accession2tg_ecotypeid ae where newt.ecotype_id is NULL and newt.accession_id=ae.accession_id) as table2, call_qc q where table2.ecotype_id=q.ecotype_id and q.qc_method_id=9 and q.call_method_id=3;

7. outer join, to get QC info, also including the ones with no QC

   select table2.*, q.call_info_id, q.QC_method_id, q.NA_rate, q.no_of_totals, q.mismatch_rate, q.no_of_mismatches, q.no_of_non_NA_pairs from (select newt.accession_id, ae.nativename, ae.stockparent, ae.ecotype_id from (select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar, call_qc q where q.call_info_id=c.id and c.method_id=3 and c.array_id=ar.id and q.qc_method_id=9 and q.mismatch_rate<0.12) order by ecotype_id, linename) as newt, at.accession2tg_ecotypeid ae where newt.ecotype_id is NULL and newt.accession_id=ae.accession_id) as table2 left outer join call_qc q on table2.ecotype_id=q.ecotype_id and q.qc_method_id=9 and q.call_method_id=3;

8. get the array id, outer join

   select t3.*, c.array_id from (select table2.*, q.call_info_id, q.QC_method_id, q.NA_rate, q.no_of_totals, q.mismatch_rate, q.no_of_mismatches, q.no_of_non_NA_pairs from (select newt.accession_id, ae.nativename, ae.stockparent, ae.ecotype_id from (select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar, call_qc q where q.call_info_id=c.id and c.method_id=3 and c.array_id=ar.id and q.qc_method_id=9 and q.mismatch_rate<0.12) order by ecotype_id, linename) as newt, at.accession2tg_ecotypeid ae where newt.ecotype_id is NULL and newt.accession_id=ae.accession_id) as table2 left outer join call_qc q on table2.ecotype_id=q.ecotype_id and q.qc_method_id=9 and q.call_method_id=3) as t3 left outer join call_info c on c.id=t3.call_info_id;

9. get array original_filename, outer join

   select t4.*, a.original_filename from (select t3.*, c.array_id from (select table2.*, q.call_info_id, q.QC_method_id, q.NA_rate, q.no_of_totals, q.mismatch_rate, q.no_of_mismatches, q.no_of_non_NA_pairs from (select newt.accession_id, ae.nativename, ae.stockparent, ae.ecotype_id from (select distinct a.accession_id, a.linename, ae.nativename, ae.stockparent, ae.ecotype_id from at.magnus_192_vs_accession a left outer join at.accession2tg_ecotypeid ae on a.accession_id=ae.accession_id and ae.ecotype_id in (select distinct ar.maternal_ecotype_id from call_info c, array_info ar, call_qc q where q.call_info_id=c.id and c.method_id=3 and c.array_id=ar.id and q.qc_method_id=9 and q.mismatch_rate<0.12) order by ecotype_id, linename) as newt, at.accession2tg_ecotypeid ae where newt.ecotype_id is NULL and newt.accession_id=ae.accession_id) as table2 left outer join call_qc q on table2.ecotype_id=q.ecotype_id and q.qc_method_id=9 and q.call_method_id=3) as t3 left outer join call_info c on c.id=t3.call_info_id) as t4 left outer join array_info a on t4.array_id=a.id order by accession_id;

4.2   2009-4-22 check missing 192 in one 250k call-method

select a.ecotype_id, nt.* from (select distinct m.*, newt.accession_id aid from (select a.accession_id, v.* from view_call v, at.accession2tg_ecotypeid a where v.call_method_id =32 and v.ecotype_id=a.ecotype_id) as newt right outer join at.magnus_192_vs_accession m on m.accession_id=newt.accession_id order by aid, accession_id) as nt, at.accession2tg_ecotypeid a where a.accession_id=nt.accession_id order by aid, accession_id;

4.3   to Provide a list of unique lines within each of the following regions

1. simply do that

   select distinct e.id, e.name, e.stockparent, e.nativename, s.name as site_name, c.abbr as country_abbr, q.call_info_id, q.mismatch_rate from stock.ecotype e, stock.site s, stock.address a, stock.country c, call_qc q where q.ecotype_id=e.id and e.siteid=s.id and s.addressid=a.id and a.countryid=c.id and q.qc_method_id=9 and q.mismatch_rate<0.12 and q.call_method_id=3 order by country_abbr;

2. remove call_info_id, mismatch_rate

   select distinct e.id as ecotype_id, e.name, e.stockparent, e.nativename, s.name as site_name, c.abbr as country_abbr from stock.ecotype e, stock.site s, stock.address a, stock.country c, call_qc q where q.ecotype_id=e.id and e.siteid=s.id and s.addressid=a.id and a.countryid=c.id and q.qc_method_id=9 and q.mismatch_rate<0.12 and q.call_method_id=3 order by country_abbr;

3. get the count for each country

   select country_abbr, count(ecotype_id) as no_of_ecotypes from (select distinct e.id as ecotype_id, e.name, e.stockparent, e.nativename, s.name as site_name, c.abbr as country_abbr from stock.ecotype e, stock.site s, stock.address a, stock.country c, call_qc q where q.ecotype_id=e.id and e.siteid=s.id and s.addressid=a.id and a.countryid=c.id and q.qc_method_id=9 and q.mismatch_rate<0.12 and q.call_method_id=3 order by country_abbr) t1 group by country_abbr order by no_of_ecotypes;

5.1   Deletion 3

It is discovered by looking at 2010 alignments as in Table 6. Despite the start-stop described in Table 8, there are 5 positions right ahead of this deletion and 5 positions right after this deletion, harboring a different haplotype compared to non-deletion haplotype.

5.2   2008-08-01 Kruskal-Wallis 250k vs 471 FRI SNPs from alignment 1843

In table at.alignment, id=1843 is around FRI. It's sequenced for 193 accessions. Pull polymorphic positions out of it as SNPs. deletion is also counted as polymorphic. Do a simple KW on it against phenotype LD. compare it with 250k (~170 accession) KW result.

Figure 2. Kruskal-Wallis 250k vs 471 FRI SNPs from alignment 1843. The gene under the 183kb-on-the-right peak, (highest in the 250k margarita result) is AT4G01040, "glycosyl hydrolase family 18 protein".

So the FRI deletions didn't pop out as strong as i would expect but considering the dependency of two deletions (Col has deletion 1 seq but not deletion 2 seq, FRI is dead still), it's sort of ok. maybe we can code two deletions as one SNP, or we just run margarita or some other haplotype-based algorithm ...

The peak within deletion 1 (around Position 268800) is ~200bp upstream of FRI's starting codon, might be some key regulatory site. The peak @268809 in deletion 1 has one crucial difference with other SNPs (such as @268785) to make it stand out. It has a base call for ecotype 8240 while SNP @268785 regard it as deletion. Due to the fact ecotype 8240's LD(long-day vernalization flowering time) is 200 (late flowering), this partial deletion 1 doesn't affect FRI function.

5.3   2008-08-04 create version 4 SNPs for 3 FRI deletions

A new alignment with same info as alignment id=1843. Pick 3 distinctive SNPs representing each deletion based on a SNP matrix figure (by DrawSNPMatrix.py).

Steps:

1. table alignment. insert a new alignment with same (chromosome, start, end) as 1843

2. table locus. create 3 new loci. as (chromosome,position,offset): (4,268809,0), (4,269962,8), (4,270712,0)

3. table allele. create 2 alleles for each locus. based on same alleles at those 3 loci

4. table genotype. copy genotype info over from 6 alignment=1843 alleles:

   insert into fri_deletion_alignment_1843 (allele, accession) select 262564, accession from genotype where allele=172949;
   insert into fri_deletion_alignment_1843 (allele, accession) select 262565, accession from genotype where allele=172950;
   insert into fri_deletion_alignment_1843 (allele, accession) select 262566, accession from genotype where allele=173485;
   insert into fri_deletion_alignment_1843 (allele, accession) select 262567, accession from genotype where allele=173484;
   insert into genotype (allele, accession, quality) select 262568, accession, quality from genotype where allele=173566;
   insert into genotype (allele, accession, quality) select 262569, accession, quality from genotype where allele=173567;
   

select f1.* , f2.* from fri_deletion_alignment_1843 f1, allele a1, fri_deletion_Shindo2005_Magnus_data_csv f2,
allele a2 where f1.accession=f2.accession and f1.allele=a1.id and f2.allele=a2.id and a1.locus=a2.locus
and f1.allele!=f2.allele order by f2.allele, f2.accession;

5.4   2008-08-05 epistasis among 3 FRI deletions

Figure 2. a primitive epistasis study among the 3 deletions in FRI (represented by '4_268809_0', '4_269962_8' and '4_270712_0'). using boolean OR to merge deletions. data is 471 SNPs discovered from alignment 1843 for 193 accessions.

Panel b,c,d are not exactly like Hagenblad2004 Figure 4, in which ecotypes are removed based on their phenotype and genotype. It has similar spirit which is that multiple genetic defects contribute to early flowering. It echoes a sense of conditional probability in bayesian networks.

Panel e,f proves the power of boolean thinking. besides, it doesn't suffer the sample size shrinking problem in conditional probability.

This figure also proves that deletion 3 in the 3 UTR region has similar effect although only 4 ecotypes have it.

5.5   2008-09-26 merge FRI deletion data with other 250k data for next-step QC/imputation

• Output two deletions with no offset (and other SNPs) by setting version=4 in variation/src/Output2010InCertainSNPs.py.:

  ./src/Output2010InCertainSNPs.py -o data/2010/data_2010_ecotype_id_y0002_version4_offset_8.tsv -r -y0002 -v 4 -f 0
  

• Remove all NA columns (only 3 deletions in version 4) by calling:

  ./src/FilterStrainSNPMatrix.py -i data/2010/data_2010_ecotype_id_y0002_version4.tsv -o 
  data/2010/data_2010_ecotype_id_y0002_version4_n1c1d110.tsv -n 1 -c 1 -d 110 -r
  

• Output the other deletion by setting version=4 and offset=8 in Output2010InCertainSNPs.py.:

  ./src/Output2010InCertainSNPs.py -o data/2010/data_2010_ecotype_id_y0002_version4_offset_8.tsv -r -y0002 -v 4 -f 8
  

• Remove all NA columns:

  ./src/FilterStrainSNPMatrix.py -i data/2010/data_2010_ecotype_id_y0002_version4_offset_8.tsv -o
  data/2010/data_2010_ecotype_id_y0002_version4_offset8_n1c1d110.tsv -n 1 -c 1 -d 110 -r
  

• identify the column with SNP_id=`4_269962` (useless because this dataset has different number of strains):

  awk '$1=="ecotype_id"{for (i=1;i<=NF;i++)  if($i=="4_269962"){print i}}' data/2010/data_2010_ecotype_id_y0002_version4_offset_8.tsv
  cut -f 176 data/2010/data_2010_ecotype_id_y0002_version4_offset_8.tsv > data/2010/4_269962_offset8.tsv
  

• merge 3 deletions all together:

  ../pymodule/TwoSNPData.py -i data/2010/data_2010_ecotype_id_y0002_version4_n1c1d110.tsv -j
  data/2010/data_2010_ecotype_id_y0002_version4_offset8_n1c1d110.tsv -o data/2010/4_269962_offset8.tsv -c1 -w1
  

• mask the deletion integer, -1, as a nucleotide code int, 2, to avoid being dropped when merging.

• 2009-1-6 add 3 FRI-deletion SNPs into 250k data matrix. Note: TwoSNPData.py removes the 2nd-column array id and replaces them with ecotype id.:

  run `pymodule/TwoSNPData.py`
  
  ~/script/pymodule/TwoSNPData.py -i ~/panfs/NPUTE_data/input/250k_l3_y.85_uniq_ecotype_20090119.tsv 
  -j ./banyan_home/script/variation/data/2010/4_269962_offset8.tsv -o 
  ~/panfs/NPUTE_data/input/250k_l3_y.85_uniq_ecotype_20090119_3_FRI_del_no_array_id.tsv  
  -c1 -w0 -m1
  

• 2009-1-27 not necessary anymore. -m1 in TwoSNPData.py carry out this role. 2009-1-6 run recoverArrayID2ndCol() from variation/src/misc.py to restore array id in the merged 250k-FRI-3-del file. WATCH: the 1st column (strain id) of input_fname has to be unique, otherwise a random one among duplicates would get a possibly wrong array id assigned.:

  input_fname =
    os.path.expanduser('~/panfs/NPUTE_data/input/250k_l3_y.85_uniq_ecotype_20080919_3_FRI_del_no_array_id.tsv')
  data_with_array_id_fname =
    os.path.expanduser('~/panfs/NPUTE_data/input/250k_l3_y.85_uniq_ecotype_20080919.tsv')
  output_fname = os.path.expanduser('~/panfs/NPUTE_data/input/250k_l3_y.85_uniq_ecotype_20080919_3_FRI_del.tsv')
  recoverArrayID2ndCol(input_fname, data_with_array_id_fname, output_fname)
  

5.6   2008-11-24 add 3 FRI-deletion SNPs into snps_context of stock_250k

They were already added into table snps on 2008-08-17 by me.

information in table snps is

| id     | name                | chromosome | position | end_position | allele1 | allele2
+--------+---------------------+------------+----------+--------------+---------+---------
| 250924 | 4_268809_0_FRI_del1 |          4 |   268809 |         NULL | A       | -      
| 250925 | 4_269962_8_FRI_del2 |          4 |   269962 |         NULL | T       | -      
| 250926 | 4_270712_0_FRI_del3 |          4 |   270712 |         NULL | A       | -     

So add them into table snps_context as touch on FRI (gene_id=828044):

snps_id    disp_pos   gene_id  gene_strand  left_or_right  disp_pos_comment
250924        64      828044  +       touch   touch
250925        936     828044  +       touch   touch
250926        1686    828044  +       touch   touch

The 1st deletion covers both promoter and initial coding region. disp_pos is calculated as gene_start_pos-deletion_stop_pos.

Found out the 3 FRI-deletion SNPs are already in table snps_context on 2008-08-19. But they only appear in papaya's db, not banyan's. what happened??

5.7   2009-9-19 Association of FLC expr (id=43) phenotype with FRI alleles as cofactor

Figure. The top two panels are LM and EMMA without any cofactor. Every two-panel beneath is one pair of comparison. The title, such as 'LM_cof_Col 6.76' , means linear model, with FRI_Col Allele as cofactor, -log(pvalue)=6.76 at that line. "cof" in the title stands for cofactor. Scores below 0.5 are omitted. NGA4 is the only gene classified as "Flower development" in GeneOntology while not in our flowering time list.

Nothing really changes (except scale). Everything we said in the paper is still valid. The increase of top pvalue in cofactor scenario is still evident. The two peaks, one around FRI, one between GA1 and DFL2, switch the top position from Col-cofactor to Ler-cofactor, which really gives us some insight into the underlying haplotype sharing (1st peak in LD with FRI-Ler allele, 2nd peak in LD with FRI-Col allele).