- Thu Jan 26 10:40:31 EST 2017
- Lucas Boatwright
import os
import pandas as pd
# Important directories
PROJ="/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare"
SCRIPTS="/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/scripts"
POLY_READS="/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/T.castellanus_trimmed_reads"
DESIGN=PROJ + "/design_files"
sorted(os.listdir(PROJ))

This step first uses bowtie to align reads and then any multiply-mapped reads and unmapped reads get passed to Last which will also attempt to uniquely align reads. Only those reads that align uniquely will be used for HSE.
align_reads_to_Crocifolius_consensed_bt_last.sbatch
align_reads_to_Lamottei_consensed_bt_last.sbatch
identify_quality.py
SeparateLast.py
parse_last_sam_v2.pl
maf-convert.py
bowtie/0.12.9
last/531
python/2.7.3
Design file:
DESIGN + trago_hybrids_with_reps.txt
Refs:
PROJ + Lamottei_post_compression_contigs.renamed.fasta
PROJ + Crocifolius_post_compression_contigs.renamed.fasta
Reads:
PROJ + /original_data/
PROJ + /alignment_output/{SP}_{REP}_uniq_2_{REF}.sam # Ex. Lamottei_1_uniq_2_Lam.sam
--Memory reduced from 160G on bigmem to 5G on HPG2. I'm not sure why the memory requirements were so large.
--After running, there is an intermediate bash script to run:
/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/filter_sam_files/filtered_parent_sam_files/rename_filtered_files.bash
In the previous step, we wanted to align to the entire transcriptome to prevent reads from mapping greedily to spurious locations in the transcriptome. Here, we isolate reads that map within our Common Orthologous REgions (COREs) so that we can make comparisons between reference mappings.
filter_sam_files_parents_2_parents_v3.sbatch
PROJ + /scripts/sam-filter-by-bed.pl
perl/5.16.0
SAM:
PROJ + /alignment_output/{SP}_{REP}_uniq_2_{REF}.sam
Bed:
PROJ + /filter_sam_files/{SP}_overlaps_WRT_orthologs.renamed.bed
PROJ + /filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_uniq_2_Lam_filter_for_Cast.sam
PROJ + /filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_uniq_2_Croc_filter_for_Cast.sam
# Ex. LAM_1_uniq_2_Lam.sam -- changed due to 32 char max length
--Recommend running directly on bbshiva as it only takes likes 3 seconds to run
Here, we import the BED files into SAS so that the BED file coordinates may be edited for use in the SAM compare step.
submit_sas_programs_4_sam_compare.sbatch
import_bed_fix_coord_for_sam_compare_v2.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
SAS scripts:
/scripts/import_bed_fix_coord_for_sam_compare_v2.sas
CSV files:
Croc-lam_overlaps_WRT_orthologs.renamed.csv
Lam-croc_overlaps_WRT_orthologs.renamed.csv
SAS output:
PROJ + /sas_data/output
PROJ + /references/Lam_croc_bed_for_sam_compare.bed
PROJ + /references/Croc_lam_bed_for_sam_compare.bed
Since the two references represent sequences from different species, we have to rename putatively orthologous pairs so that SAM compare is comparing the correct sequence pairs.
submit_sas_program_add_commonID_parent.sbatch
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
SAS scripts:
add_commonID_to_parent_sam_files_HPC.sas
/filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_uniq_2_Lam.sam
/filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_uniq_2_Croc.sam
SAS output:
PROJ + /sas_data/output
/filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_unq_2_Lam_commonID.sam
/filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_unq_2_Croc_commonID.sam
This step performs SAM compare to obtain counts for each reference.
sam_compare_parental_reads.sbatch
sam_compare.py
python/2.7.3
/filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_unq_2_Lam_commonID.sam
/filter_sam_files/filtered_parent_sam_files/{SP}_{REP}_unq_2_Croc_commonID.sam
/ase_counts_{SP}_{REP}_2_Lam_Croc.csv
/ase_totals_{SP}_{REP}_2_Lam_Croc.csv
-- SAS scripts located at PROJ + /sas_downstream_scripts/
-- I actually ran this on bbshiva...takes like 3s
This step is self-explanatory
submit_import_ase_counts_parents.bash
import_ase_counts_parents.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
/ase_counts_{SP}_{REP}_2_Lam_Croc.csv
PROJ + /sas_data/output
--edited for parents and to correct for change above
This step takes a matrix of counts:
total_ ase_ Lam_ Croc_
Obs commonID rep count count ase count count hybrid
1 Lamottei_100001|Crocifolius_2415 1 40 22 0.95455 21 1 Lam
2 Lamottei_100001|Crocifolius_2415 2 26 17 1.00000 17 0 Lam
3 Lamottei_100001|Crocifolius_2415 3 41 26 1.00000 26 0 Lam
4 Lamottei_100001|Crocifolius_2415 4 51 31 1.00000 31 0 Lam
5 Lamottei_100001|Crocifolius_2415 5 46 29 1.00000 29 0 Lam
6 Lamottei_100001|Crocifolius_2415 6 45 25 0.00000 0 25 Cro
7 Lamottei_100001|Crocifolius_2415 7 53 30 0.00000 0 30 Cro
8 Lamottei_100001|Crocifolius_2415 8 62 37 0.00000 0 37 Cro
9 Lamottei_100003|Crocifolius_5465 1 70 70 1.00000 70 0 Lam
10 Lamottei_100003|Crocifolius_5465 2 44 43 1.00000 43 0 Lam
And places the reps into individual columns (side-by-side).
submit_bayesian_make_sbys_reps_parents.bash
bayesian_make_sbys_reps_parents.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
SAS data-
trago.ase_4_bayes_&reads._reads_lam_croc
/sas_data/output/
This step flags each feature (putative orthologous pair) for analysis (1) or not (0). Each feature MUST have at least one read and all reps should be present. This check for presence is important when running different rep counts.
submit_bayesian_flag_analyze_hybrid.bash
bayesian_flag_analyze_parents.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
SAS data-
trago.ase_4_bayes_cast_reads_lam_croc
SAS data-
trago.bayes_flag_cast_reads_lam_croc
Self-explanatory -- export those data which were flagged for analysis (1) so that we can run the PG analysis (Empirical Bayesian Machine).
submit_bayesian_export_flagged_dataset_parents.bash
bayesian_export_flagged_dataset_parents.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
ase_bayes_par_sbys_with_flag
PROJ + /empirical_bayesian_hybrid_input/ase_bayes_par_lam_croc_flag.csv
-- Had to comment out some logging commands that had a dependency I was missing -- not an issue
This step is performed to speed up the analysis by separating the process into reasonably sized chunks (they take about 4 minutes or less to run)
split_emp_bayesian_ase_table.sh
splitTable.py
python/2.7.6
/empirical_bayesian_parents_input/ase_bayes_croc_lam_croc_flag.csv
/empirical_bayesian_parents_input/ase_bayes_lam_lam_croc_flag.csv
/empirical_bayesian_parents_input/output/split_{SP}/split_1.csv .. 500
Execute the Bayesian PG machine to identify expression bias between the orthologous pairs. For the parents, parental reads were mapped back to both references to identify loci that exhibit a mapping bias toward the wrong parental reference.
run_emp_bayesian_machine_parents_PG5.sbatch
AI_poissongamma_functions.R
PG_model_empirical_q456.r
Subroutines_model2_experimental.R
# Additional scripts as PG only works for 3 reps
PG_model_empirical_q456_5reps.r
R/3.0.2
/empirical_bayesian_parents_input/output/split_{SP}/split_1.csv .. 500
empirical_bayesian_parents_output/split_{SP}/split_1.csv .. 500
-- Had to comment out some logging commands that had a dependency I was missing -- not an issue
After each of the split runs of the Emp_Bayes_Machine finishes, all of the output is subsequently merged back into one file.
combine_emp_bayesian_ase_table_parents.sh
catTable.py
python/2.7.6
empirical_bayesian_parents_output/split_{SP}/split_1.csv .. 500
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents.csv
Self-explanatory: Import the EBM output into SAS.
execute_import_bayesian_results.bash
import_bayesian_results.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents.csv
/sas_data/
Flag those putative orthologous pairs that exhibit mapping bias. Flag significant if the
credible interval [similar to 95% Highest posterior density interval (HPDI)] does not overlap 0.5
for each prior (q4,q5,q6).
Then, flag significant for all three if q4 = q5 = q6 = 1
execute_bayesian_flag_sig_results_parents.bash
bayesian_flag_sig_results_parents_pdf.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents.csv
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents_flag_sig.csv
from IPython.display import Image
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/bayes_flag_sig_Lam_for_Lam_Croc.png')
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/bayes_flag_sig_Croc_for_Lam_Croc.png')

- We have 4 independent T. castellanus polyploids with 3 bioreps each
- We really only have 1 population for both T. lamottei and T. crocifolius
- We need to assess HSE in each T. castellanus population separately
- NOTE: We have parents and polyploid P-B + P-II = IIb --> we should assess HSE specially for these
Image(filename='T.castellanus_parentage.jpg')
# NOTES:
# Graph not entirely accurate
# T.cast IIa has an additive karyotype
# T.cast IIb has a reciprocal translocation --> this is what we have parent data for
# How samples compare to parentage above
samples = pd.read_csv("T.croc_T.lam_samples.csv", index_col="Sample ID")
samples
pre-emptively split scripts since the perl portion takes bigmem to properly execute
this took a ton of memory -- especially for alignment to Croc
align_Castellanus_reads_to_Lamottei_consensed_bt_last.sbatch
align_Castellanus_reads_to_Crocifolius_consensed_bt_last.sbatch
align_Castellanus_reads_to_Lamottei_perl_portion.sbatch
align_Castellanus_reads_to_Crocifolius_perl_portion.sbatch
identify_quality.py
SeparateLast.py
parse_last_sam_v2.pl
maf-convert.py
bowtie/0.12.9
last/531
python/2.7.3
Design file:
DESIGN + trago_cast_with_reps.txt
Refs:
PROJ + {SP}_post_compression_contigs.renamed.bowtie.ref
PROJ + {SP}_post_compression_contigs.renamed.last531.ref
Reads:
POLY_READS + /{SP}{REP}/*.fq.gz
PROJ + /alignment_output/{SP}_{REP}_uniq_2_{REF}.sam # Ex. Lamottei_1_uniq_2_Lam.sam
filter_sam_files_cast_2_parents_v3.sbatch
PROJ + /scripts/sam-filter-by-bed.pl
perl/5.16.0
SAM:
PROJ + /alignment_output/{SP}_{REP}_uniq_2_{REF}.sam
Bed:
PROJ + /filter_sam_files/{SP}_overlaps_WRT_orthologs.renamed.bed
PROJ + /filter_sam_files/filtered_parent_sam_files/Cast_{REP}_uniq_2_Lam_filter_for_Cast.sam
PROJ + /filter_sam_files/filtered_parent_sam_files/Cast_{REP}_uniq_2_Croc_filter_for_Cast.sam
# Ex. LAM_1_uniq_2_Lam.sam -- changed due to 32 char max length?
python rename_cast_files.py /ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/filter_sam_files/filtered_parent_sam_files
rename_cast_files.py
Cast_13-3-5_uniq_2_Lam_filter_for_Cast.sam
Cast_13_3_5_uniq_2_Lam_filter_for_Cast.sam
submit_add_commonID_hybrid.bash
sas/9.4
SAS scripts:
add_commonID_to_hybrid_sam_files_HPC.sas
/filter_sam_files/filtered_parent_sam_files/Cast_{REP}_uniq_2_Lam.sam
/filter_sam_files/filtered_parent_sam_files/Cast_{REP}_uniq_2_Croc.sam
SAS output:
PROJ + /sas_data/output
/filter_sam_files/filtered_parent_sam_files/Cast_{REP}_unq_2_Lam_commonID.sam
/filter_sam_files/filtered_parent_sam_files/Cast_{REP}_unq_2_Croc_commonID.sam
sam_compare_hybrid_reads.sbatch
sam_compare.py
python/2.7.3
/filter_sam_files/filtered_parent_sam_files/Cast_{REP}_unq_2_Lam_commonID.sam
/filter_sam_files/filtered_parent_sam_files/Cast_{REP}_unq_2_Croc_commonID.sam
/ase_counts_Cast_{REP}_2_Lam_Croc.csv
/ase_totals_Cast_{REP}_2_Lam_Croc.csv
This step and all following steps are located at PROJ/sas_downstream_scripts/
Cast CSV files were copied and renamed due to SAS macro issues with underscores and hyphens...
submit_import_ase_counts_hybrid.bash
import_ase_counts_hybrids_renamed.sas
sas/9.4
/ase_counts_{SP}_{REP}_2_Lam_Croc.csv (ex. ase_counts_Cast_242_2_Lam_Croc.csv)
PROJ + /sas_data/output
-- Had to comment out some logging commands that had a dependency I was missing -- not an issue
split_emp_bayesian_ase_table_hybrids.bash
splitTable.py
python/2.7.6
/empirical_bayesian_parents_input/ase_bayes_cast2_lam_croc_flag.csv
.
.
/empirical_bayesian_parents_input/ase_bayes_cast31_lam_croc_flag.csv
/empirical_bayesian_parents_input/output/split_{SP}/split_1.csv .. 500
Execute the Bayesian PG machine to identify expression bias between the orthologous pairs. For the hybrid, reads were mapped back to both references to identify loci that exhibit an expression bias toward each parental reference. This was done using all loci and is later filtered (keeping loci that do not exhibit mapping bias -- derived from the parent PG pipe).
run_emp_bayesian_machine_hybrids.sbatch
AI_poissongamma_functions.R
PG_model_empirical_q456.r
Subroutines_model2_experimental.R
R/3.0.2
/empirical_bayesian_hybrid_input/output/split_{SP}/split_1.csv .. 500
empirical_bayesian_hybrid_output/split_{SP}/split_1.csv .. 500
-- Had to comment out some logging commands that had a dependency I was missing -- not an issue
combine_emp_bayesian_ase_table_hybrid.bash
catTable.py
python/2.7.6
empirical_bayesian_hybrid_output/split_{SP}/split_1.csv .. 500
empirical_bayesian_hybrid_output/PG_emp_bayesian_results_{SP}_hybrid.csv
execute_import_bayesian_flag_results_hybrid.bash
import_bayesian_results_hybrid.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
empirical_bayesian_hybrid_output/PG_emp_bayesian_results_{SP}_hybrid.csv
/sas_data/
execute_bayesian_flag_sig_results_hybrid.bash
bayesian_flag_sig_results_hybrid_pdf.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents.csv
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents_flag_sig.csv
from IPython.display import Image
print("Plot of mean theta for all POPs (putative orthologous pairs) -- whether significantly biased or not")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot.png')
from IPython.display import Image
print("Plot of mean theta for all POPs (putative orthologous pairs) -- whether significantly biased or not")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot1.png')
from IPython.display import Image
print("Plot of mean theta for all POPs (putative orthologous pairs) -- whether significantly biased or not")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot2.png')
from IPython.display import Image
print("Plot of mean theta for all POPs (putative orthologous pairs) -- whether significantly biased or not")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot3.png')
from IPython.display import Image
print("Plot of mean theta for POPs (putative orthologous pairs) -- significantly biased using all priors")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot4.png')
from IPython.display import Image
print("Plot of mean theta for POPs (putative orthologous pairs) -- significantly biased using all priors")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot5.png')
from IPython.display import Image
print("Plot of mean theta for POPs (putative orthologous pairs) -- significantly biased using all priors")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot6.png')
from IPython.display import Image
print("Plot of mean theta for POPs (putative orthologous pairs) -- significantly biased using all priors")
Image(filename='/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/reports/SIG3_hybrid_plots/SGPlot7.png')

-- SAS scripts located at PROJ + /sas_downstream_scripts/
-- I actually ran this on bbshiva...takes like 3s
submit_import_ase_counts_parents_Lam3.bash
import_ase_counts_parents_Lam3.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
/ase_counts_{SP}_{REP}_2_Lam_Croc.csv
PROJ + /sas_data/output
--edited for parents and to correct for change above
submit_bayesian_make_sbys_reps_parents_L3.bash
bayesian_make_sbys_reps_parents_L3.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
SAS data-
trago.ase_4_bayes_&reads._reads_lc_l3
/sas_data/output/ase_bayes_&reads._rds_lc_l3_sbys
submit_bayesian_flag_analyze_parents_L3.bash
bayesian_flag_analyze_parents_L3.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
SAS data-
trago.ase_4_bayes_cast_reads_lc_l3
SAS data-
trago.bayes_flag_cast_reads_lam_croc
submit_bayesian_export_flagged_dataset_parents_L3.bash
bayesian_export_flagged_dataset_parents_L3.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
Bayes_flag_&reads._reads_lc_l3
Ase_bayes_&reads._rds_lc_l3_sbys
PROJ + /empirical_bayesian_hybrid_input/ase_bayes_{SP}_lc_l3_flag.csv
-- Had to comment out some logging commands that had a dependency I was missing -- not an issue
--Move back to 'scripts' directory
split_emp_bayesian_ase_table_L3.bash
splitTable.py
python/2.7.6
/empirical_bayesian_parents_input/ase_bayes_croc_lc_l3_flag.csv
/empirical_bayesian_parents_input/ase_bayes_lam_lc_l3_flag.csv
/empirical_bayesian_parents_input/output/split_{SP}/split_1.csv .. 500
Move to PROJ/scripts/emp_bayesian/
run_emp_bayesian_machine_parents_L3.sbatch
AI_poissongamma_functions.R
PG_model_empirical_q456.r
Subroutines_model2_experimental.R
R/3.0.2
/empirical_bayesian_parents_input/output/split_{SP}_l3/split_1.csv .. 500
empirical_bayesian_parents_output/split_{SP}/split_1_l3.csv .. 500
-- Had to comment out some logging commands that had a dependency I was missing -- not an issue
combine_emp_bayesian_ase_table_parents_L3.bash
catTable.py
python/2.7.6
empirical_bayesian_parents_output/split_{SP}_l3/split_1.csv .. 500
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents_l3.csv
In sas_downstream directory
execute_import_bayesian_flag_results_L3.bash
import_bayesian_results_L3.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents.csv
/sas_data/
execute_bayesian_flag_sig_results_parents_L3.bash
bayesian_flag_sig_results_parents_pdf_L3.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents_l3.csv
empirical_bayesian_parents_output/PG_emp_bayesian_results_{SP}_parents_flag_sig_l3.csv
from IPython.display import Image
Image("./Lam_l3.png")
Image("Croc_l3.png")

Export the flagged EBM results into a CSV for filtering.
execute_output_ase_bayesian_flagged_CSVs.bash
output_ase_bayesian_flagged_CSVs.sas
sas/9.4 (originally 9.3 but due to licensing, this must be updated accordingly)
/sas_data/trago.bayes_flag_sig_Cast2
/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/bayes_flag_sig_Cast2.csv
This step identifies parental sequences that biasedly mapped reads from the 'wrong' parent and uses the unbiased set of loci to filter the HSE hybrid data.
None, but there are several functions and commands below.
python libraries:
pandas
numpy
matplotlib
matplotlib_venn
python/2.7.6
/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/bayes_flag_sig_{SP}_for_LC_L3.csv
Figures and filtered hybrid HSE CSV
def filter_csv(file_name):
"""Given flagged bayesian output in CSV format, isolate unbiased loci.
This is only meant for the parents, as the parent data can biasedly map to its own reference or exhibit no bias."""
import pandas as pd
import numpy as np
df = pd.read_csv(file_name, index_col = "commonID")
cols = list(df.columns)
sig_col = cols[-1]
reads=sig_col.split("_")[2]
unbiased_loci = []
print("\tDF pre-filter: " + str(len(df.index)))
if reads == "Lam":
# identify loci not biased toward Croc
not_biased = df[(df[sig_col] == 0) & (df["q5_mean_theta"] <= 0.75)]
print("\t\tNot biased: " + str(len(not_biased.index)))
df = pd.concat([df[(df["q5_mean_theta"] <= 0.5) & (df[sig_col] == 1)],not_biased])
elif reads == "Croc":
# identify loci not biased toward Lam
not_biased = df[(df[sig_col] == 0) & (df["q5_mean_theta"] >= 0.25)]
print("\t\tNot biased: " + str(len(not_biased.index)))
df = pd.concat([df[(df["q5_mean_theta"] >= 0.5) & (df[sig_col] == 1)], not_biased])
print("\tDF post-filter: " + str(len(df.index)))
return df
def inline_hist(df, file_name):
"""Inline plot histogram of q5 mean theta"""
plt.figure()
plt.hist(df["q5_mean_theta"])
plt.title(file_name + " q0.5 mean theta")
plt.xlabel("Theta")
plt.ylabel("Frequency")
plt.show()
def csv_to_df(file_name):
"""Given flagged bayesian output in CSV format, convert to standard dataframe.
This is only meant for the hybrid, as the hybrid data may biasedly express one ortholog over another."""
import pandas as pd
import numpy as np
df = pd.read_csv(file_name, index_col = "commonID")
cols = list(df.columns)
cols[-1]="sig_all_three"
df.columns = cols
print("Length of dataframe pre-filter: " + str(len(df["sig_all_three"])))
return df
def hist_unfiltered_filtered(df, unbiased_loci):
"""Plot a histogram of both the unfiltered and filtered q5 mean theta to visualize filtering"""
import numpy as np
import matplotlib.pyplot as plt
plt.figure()
unbiased_df = df[df.index.isin(unbiased_loci)]
print("\tLoci exhibiting HSE pre-filter: " + str(np.sum(df["sig_all_three"])))
n, bins, patches = plt.hist(df["q5_mean_theta"], bins=50, histtype='step', facecolor='r', label='pre-filter')
print("Total loci post-filter: " + str(len(unbiased_df["sig_all_three"])))
print("\tLoci exhibiting HSE post-filter: " + str(np.sum(unbiased_df["sig_all_three"])))
n, bins, patches = plt.hist(unbiased_df["q5_mean_theta"], bins=50, histtype='step', facecolor='b', label='post-filter')
plt.xlabel('Mean Theta')
plt.ylabel('Frequency')
plt.legend()
plt.title('Distribution of Mean Theta for T. castellanus')
plt.rcParams['figure.figsize'] = (8.0, 8.0)
plt.show()
return unbiased_df
import os
file_path = "/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/"
files = [i for i in os.listdir(file_path) if "LC" in i]
df_dict = {}
for i in files:
print(i) # name of the file
df_dict[i]= filter_csv(file_path + i) # Not biased just means not biased toward wrong parent (bias toward self is okay)
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib_venn import venn2_unweighted, venn2_circles
# plot the commonID overlap after filtering for biased loci
# The overlap represents the set of loci for downstream hybrid analysis
set_Croc = set(df_dict["bayes_flag_sig_Croc_for_LC.csv"].index)
set_Lam = set(df_dict["bayes_flag_sig_Lam_for_LC.csv"].index)
plt.figure()
venn2_unweighted(subsets=[set_Lam, set_Croc], set_labels=('Lam','Croc'))
plt.show()
# print histograms of filtered mean theta values to assess distribution -- hard cutoffs at 0.75 and 0.25
# These plots indicate that sequences are allowed to map more reads from their corresponding parent
for key, value in df_dict.items():
if (("Lam" in key) and ("LC." in key)):
lam_lc = df_dict[key][df_dict[key].index.isin(set_Lam)]
inline_hist(lam_lc, key)
elif (("Croc" in key) and ("LC." in key)):
croc_lc = df_dict[key][df_dict[key].index.isin(set_Croc)]
inline_hist(croc_lc, key)
import os
file_path = "/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/"
# get list of Cast (hybrid) files
files = [i for i in os.listdir(file_path) if (("Cast" in i) and ("putative" not in i))]
# pull out unbiased loci
unbiased_LC_loci = set_Croc.intersection(set_Lam)
with open("unbiased_LC_loci.txt",'w') as f:
for locus in list(unbiased_LC_loci):
f.write(locus + "\n")
cast_df_dict = {}
unbiased_cast_dfs = {}
for file in files:
# Read data into dataframe
df = csv_to_df(file_path + file)
# generate histogram of Cast data pre- and post-filter
unbiased_cast_dfs[file] = hist_unfiltered_filtered(df, unbiased_LC_loci)
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib_venn import venn2_unweighted, venn2_circles
# plot the commonID overlap after filtering for biased loci
set_Croc = set(df_dict["bayes_flag_sig_Croc_for_LC_L3.csv"].index)
set_Lam = set(df_dict["bayes_flag_sig_Lam_for_LC_L3.csv"].index)
plt.figure()
venn2_unweighted(subsets=[set_Lam, set_Croc], set_labels=('Lam','Croc'))
plt.show()
# print histograms of filtered mean theta values to assess distribution -- hard cutoffs at 0.75 and 0.25
# These plots indicate that sequences are allowed to map more reads from their corresponding parent
for key, value in df_dict.items():
if (("Lam" in key) and ("LC_L3" in key)):
lam_lc = df_dict[key][df_dict[key].index.isin(set_Lam)]
inline_hist(lam_lc, key)
elif (("Croc" in key) and ("LC_L3" in key)):
croc_lc = df_dict[key][df_dict[key].index.isin(set_Croc)]
inline_hist(croc_lc, key)
import os
file_path = "/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/"
files = [i for i in os.listdir(file_path) if (("Cast" in i) and ("putative" not in i))]
unbiased_L3_loci = set_Croc.intersection(set_Lam)
with open("unbiased_L3_loci.txt",'w') as f:
for locus in list(unbiased_L3_loci):
f.write(locus + "\n")
cast_df_dict = {}
unbiased_cast_dfs_l3 = {}
for file in files:
df = csv_to_df(file_path + file)
unbiased_cast_dfs_l3[file] = hist_unfiltered_filtered(df, unbiased_L3_loci)
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib_venn import venn2_unweighted, venn2_circles
for key, value in unbiased_cast_dfs.items():
print(key)
lc = set(value.index)
l3 = set(unbiased_cast_dfs_l3[key].index)
# plot the commonID overlap for LC and L3
plt.figure()
venn2_unweighted(subsets=[lc, l3], set_labels=('LC','L3'))
plt.show()
import os
os.chdir("./unbiased_overlap/")
for key in unbiased_cast_dfs_l3.keys():
output = key.split(".csv")[0] + "_unbiased_L3.csv"
unbiased_cast_dfs_l3[key].to_csv(output)
!Rscript UpSetR_unbiased_L3.R
from IPython.display import Image
Image("upset_unbiased_plot.png")
## Overlap in Cast biased -- note: bias direction not assessed in this preliminary check
list_of_cast_filtered_dfs = [set(df[df["sig_all_three"]==1].index) for df in unbiased_cast_dfs.values()]
common_biased = set.intersection(*list_of_cast_filtered_dfs)
print(len(common_biased))
This step identifies putative homeolog loss from the set of unbiased contigs. The set is first filtered by looking at homeolog pairs where one reference has zero reads. Then, there must be HSE toward the parent with data (i.e. the posterior may not overlap theta 0.5).
putative_homeolog_loss.py
python libraries:
pandas
python/2.7.6
unbiased references
./empirical_bayesian_hybrid_input/ase_bayes_Cast#_lam_croc_flag.csv -- for counts
./ase_results/bayes_flag_sig_Cast2.csv -- for HSE
Putatively lost homeologs
# %load /ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/scripts/putative_homeolog_loss.py
%%writefile /ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/scripts/putative_homeolog_loss.py
"""
Created on Tue May 9 13:48:20 2017
@author: Lucas Boatwright
"""
from sys import argv, stderr
import argparse
from datetime import datetime
import pandas as pd
def parse_arguments():
"""Parse arguments passed to script"""
parser = argparse.ArgumentParser(description=
"This script is designed to check for putative loss between homeologs \
\nusing counts from homeolog COREs used in the PG analysis.\n\n \
Example: python {0} -c ase_bayes_Cast2_lam_croc_flag.csv -b bayes_flag_sig_Cast2.csv".format(argv[0]),
formatter_class = argparse.RawDescriptionHelpFormatter)
requiredNamed = parser.add_argument_group('required arguments')
requiredNamed.add_argument("-c", "--COUNTS", type=str, required=True,
help="Bayes flagged CSV of read counts for Line and Tester", action="store")
requiredNamed.add_argument("-b", "--BIAS", type=str, required=True,
help="Bayes flagged CSV of loci demonstrating HSE", action="store")
return parser.parse_args()
def find_putative_loss(file_name):
"""Given file name, determine loci that are putatively lost based upon lack of expression"""
df = pd.read_csv(file_name)
output_base = file_name.split(".csv")[0]
df["Line_mean"] = df[["LINE_TOTAL_1","LINE_TOTAL_2","LINE_TOTAL_3"]
].mean(axis=1)
df["Line_std"] = df[["LINE_TOTAL_1","LINE_TOTAL_2","LINE_TOTAL_3"]
].std(axis=1)
df["Tester_mean"] = df[["TESTER_TOTAL_1","TESTER_TOTAL_2","TESTER_TOTAL_3"]
].mean(axis=1)
df["Tester_std"]=df[["TESTER_TOTAL_1","TESTER_TOTAL_2","TESTER_TOTAL_3"]
].std(axis=1)
df["Line_zero"] = ((df["Tester_mean"] > 0 ) & (df["Line_mean"] == 0 ))
df["Tester_zero"] = ((df["Line_mean"] > 0 ) & (df["Tester_mean"] == 0 ))
return df
def verify_HSE(df, hse_file):
"""Verify that each locus does not overlap 0.5 based upon 95% HDI"""
flag_sig = pd.read_csv(hse_file)
sig_true = flag_sig.columns[-1]
merged_df = pd.merge(df, flag_sig)
output_df = merged_df[["commonID","Line_zero","Tester_zero","Line_mean","Tester_mean",sig_true]]
lost = output_df[(output_df[sig_true]==1) & (output_df["Line_zero"] | output_df["Tester_zero"])]
output_base = hse_file.split(".csv")[0] + "_putative_loss.csv"
lost.to_csv(output_base, index=False)
if __name__ == "__main__":
start = datetime.now()
args = parse_arguments()
stderr.write("Executed: python {0} -c {1} -b {2}\n".format(argv[0],args.COUNTS, args.BIAS))
df = find_putative_loss(args.COUNTS)
verify_HSE(df, args.BIAS)
stop = datetime.now()
stderr.write("Runtime: {0}\n".format(str(stop - start)))
import os
os.chdir("/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/")
files = [i for i in os.listdir("./empirical_bayesian_hybrid_input/") if ".csv" in i]
files
%%bash
PROJ="/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare"
python ${PROJ}/scripts/putative_homeolog_loss.py \
-c ${PROJ}/empirical_bayesian_hybrid_input/ase_bayes_Cast2_lam_croc_flag.csv \
-b ${PROJ}/ase_results/bayes_flag_sig_Cast2.csv
python ${PROJ}/scripts/putative_homeolog_loss.py \
-c ${PROJ}/empirical_bayesian_hybrid_input/ase_bayes_Cast10_lam_croc_flag.csv \
-b ${PROJ}/ase_results/bayes_flag_sig_Cast10.csv
python ${PROJ}/scripts/putative_homeolog_loss.py \
-c ${PROJ}/empirical_bayesian_hybrid_input/ase_bayes_Cast13_lam_croc_flag.csv \
-b ${PROJ}/ase_results/bayes_flag_sig_Cast13.csv
python ${PROJ}/scripts/putative_homeolog_loss.py \
-c ${PROJ}/empirical_bayesian_hybrid_input/ase_bayes_Cast31_lam_croc_flag.csv \
-b ${PROJ}/ase_results/bayes_flag_sig_Cast31.csv
ls /ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/*putative*
!head /ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results/bayes_flag_sig_Cast10_putative_loss.csv
%%bash
cd /ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results
module add R
Rscript UpSetR_loss_L3.R
import os
os.chdir("/ufrc/barbazuk/lboat/T.lamottei_T.crocifolius/Complete_data_set/align_reads_for_samcompare/ase_results")
from IPython.display import Image
Image("upset_loss_plot.png")


In simple terms, we do mostly theoretical research in kinematics, dynamics and control. Sort of like that deals with solving very long equations.
Know More »
We are developing an endoscopy device in line with a snake robot.
We are developing an optimum model for solar tracker.


@* * Copyright 2010-2017 Artima, Inc. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. *@ @(selectedStyle: examples.StyleTraitExamples) @import controllers.Application.milestoneScaladoc
import org.scalatest._ import OptionValues._ val anOption = Some(18) // Use OptionValues to check Option's value, using assert assert(anOption.value > 9) import Matchers._ val anOption = Some(18) // Use OptionValues to check Option's value, using Matchers syntax anOption.value should be > 9