<jupyter_start><jupyter_text>Notebook Analyzes Duplication in Motif Defintions<jupyter_code>import venusar
import motif
import thresholds
import motifs
import activity
import tf_expression
import gene_expression
# to get code changes
import imp  
imp.reload(motif)<jupyter_output><empty_output><jupyter_text>Import Motif Definitions<jupyter_code>motif_f_base='../../data/HOCOMOCOv10.JASPAR_FORMAT.TF_IDS.txt'
pc = 0.1
th = 0
bp = [0.25, 0.25, 0.25, 0.25]
motif_set_base = motif.get_motifs(motif_f_base, pc, th, bp)
motif_set_base.length()
motif_dict = motif_set_base.motif_count(False)
len(motif_dict)

# drop to valid motifs only; do twice to see if any invalid
motif_dict = motif_set_base.motif_count(True)
len(motif_dict)
range(len(motif_set_base.motifs))
range(0, (len(motif_set_base.motifs) - 1))<jupyter_output><empty_output><jupyter_text>Convert to DataFrame for Analysis of Duplication<jupyter_code>import pandas

# -- building data frame from dictionary
#df = pandas.DataFrame(motif_dict) # errors about 'you must pass an index'
# ref: http://stackoverflow.com/questions/17839973/construct-pandas-dataframe-from-values-in-variables#17840195
dfCounts = pandas.DataFrame(motif_dict, index=[0])
dfCounts = pandas.melt(dfCounts) # rotate
dfCounts.rename(columns = {'variable': 'TF', 'value': 'tfCount'}, inplace=True)

# -- handling types
# ref: http://stackoverflow.com/questions/15891038/pandas-change-data-type-of-columns
dfCounts.dtypes

#pandas.to_numeric(s, errors='ignore')

dfCounts.size
dfCounts[dfCounts.tfCount > 1]
# add additional information to the dataframe
msb_names,msb_lengths = motif_set_base.motif_lengths(False)
# define data frame by columns
dfc = pandas.DataFrame(
    {
        "TF" : msb_names,
        "TFLength" : msb_lengths
    })
print(dfc)

if False:
    # define data frame by row (this is wrong; interlaced value sets)
    dfr = pandas.melt(pandas.DataFrame(
        [ msb_names,
        msb_lengths])).rename(columns = {'variable': 'TF', 'value': 'tfCount'}, inplace=True)
    dfr
dfCounts[dfCounts.tfCount > 1].TF
dfc
duplication_set = pandas.merge( dfCounts[dfCounts.tfCount > 1], dfc, how='inner', on='TF' ).sort_values('TF')
duplication_set
#duplication_set.select(['TF','TFLength']).groupby(by='TF').rank(method='min')
duplication_set.groupby(by='TF').rank(method='min')
# ref: http://stackoverflow.com/questions/23976176/ranks-within-groupby-in-pandas

dfRank = lambda x: pandas.Series(pandas.qcut(x,2,labels=False),index=x.index)
dfRank2 = lambda x: pandas.qcut(x,2,labels=False)
# this works: replacing x above with duplication_set['TFLength'] but fails when adding groupby, why? fails using apply too.
# duplication_set['ranks'] = duplication_set.groupby('TF')['TFLength'].apply(dfRank)
# duplication_set['ranks'] = duplication_set['TFLength'].apply(dfRank)
# duplication_set['ranks'] = pandas.qcut((duplication_set['TFLength']),2,labels=False)
duplication_set['ranks'] = dfRank2(duplication_set['TFLength'])

# adding rank to try to pivot multiple rows to columns but no dice
# df.pivot(columns='var', values='val')
#duplication_set[['TF','TFLength','ranks']].pivot(columns='ranks',values='TFLength') # stupidly keeps dropping TF column, why? also duplicating rows and not actually pivoting
# duplication_set.pivot_table(df,index=["TF","Ranks"]) # fails, 'grouper for TF not 1 dimensional
# hack tired of fighting odd outcome, dyplr is much better than pandas
duplication_set

pandas.qcut(duplication_set['TFLength'],2,labels=False)  # doesn't error gives ranks

# this is wrong but not clear why?
duplication_set[['TF','TFLength','ranks']].pivot(index='TF', columns='ranks',values='Lengths') # this is wrong too

duplication_set

duplication_set.dtypes

# duplication_set.pivot(index='TF', columns='ranks',values='TFLength') # this is wrong too errors: pandas pivot ValueError: Index contains duplicate entries, cannot reshape
# led to
# ref: http://stackoverflow.com/questions/28651079/pandas-unstack-problems-valueerror-index-contains-duplicate-entries-cannot-re#28652153
# e.set_index(['id', 'date', 'location'], append=True)
# not this doesn't work either and just creates problems
#duplication_set.set_index(['TF', 'ranks', 'TFLength','tfCount'], append=True)#.pivot( columns='ranks',values='TFLength')
duplication_set[duplication_set.TF == 'RFX5']<jupyter_output><empty_output><jupyter_text>Move On: Just focus on Direct Runs- Debugging the motifs.py code; - grep 'RFX5' output.motif.20170114.vcf | grep '5.4636' > temp_investigate.txt- duplication of RFX5 occurs in chr1, 145039922; chr22 25005718; and chr5 46363751 Build up Reduced Data Set Built reduced Motif File:- only include RFX5- note the motif file must have an end blank line- name: HOCOMOCOv10.JASPAR_FORMAT.TF_IDS.fpr_0p001.RFX5.txt Building Example VCF File [Code]- grep '^' ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.vcf > ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf - grep 'chr1' ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.vcf | grep '145039922' >> ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf- grep 'chr22' ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.vcf | grep '25005718' >> ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf- grep 'chr5' ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.vcf | grep '46363751' >> ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf building example vcf file for ZNF143 [Code]- grep '^' ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.vcf > ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset2.vcf- grep 'chr1' ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.vcf | grep '762601' >> ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset2.vcf Running the Replacement Code- python3 tf_expression.py -i ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf -o 1 -m ../../data/HOCOMOCOv10.JASPAR_FORMAT.TF_IDS.fpr_0p001.RFX5.txt -e ../../data/ALL_ARRAYS_NORMALIZED_MAXPROBE_LOG2_COORDS.sorted.txt -mo ../../data/HOCOMOCOv10.JASPAR_FORMAT.TF_IDS.fpr_0p001.txt.bed_reduced.RFX5.txt    + Creating gene dictionary for expression data.    + start read exp_file:../../data/ALL_ARRAYS_NORMALIZED_MAXPROBE_LOG2_COORDS.sorted.txt    + Filtering motif info for TFs that don't meet the expression threshold of 5. Found 17225 genes. Start filtering motifs.    + COMPLETE.- python3 motifs.py -i ../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf -r ../../data/genome_reference/reference_genome_hg19.fa -m ../../data/HOCOMOCOv10.JASPAR_FORMAT.TF_IDS.fpr_0p001.txt.bed_reduced.RFX5.txt -o ../../data/output.motif.20170507.vcf -fm -fp -ci ../../data/GM12878.ENCODE.ALL_TFS.bed -co ../../data/output.chip_peaks_output.20170507.bed &> ../../data/0_run_logs/20170507_motifs_run_stdout.txt Examining Logs- Shows 2 matches, 1 positive, 1 negative (ie reverse complement)-<jupyter_code># repeating motifs.py sub code set

import vcf
import motifs
import sequence
from pyfaidx import Fasta

force_ref_match = False

file_motif='../../data/HOCOMOCOv10.JASPAR_FORMAT.TF_IDS.fpr_0p001.txt.bed_reduced.RFX5.txt'
pc = 0.1
th = 0
bp = [0.25, 0.25, 0.25, 0.25]
ws = 50
motif_set = motif.get_motifs(file_motif, pc, th, bp)
wing_l = max(motif_set.max_positions, ws)

file_reference_genome='../../data/genome_reference/reference_genome_hg19.fa'
fa_ind = Fasta(file_reference_genome)    # XXX: need to check, if present skip


file_input = '../../data/FLDL_CCCB_RARE_VARIANTS.MERGED.RNA_DP10.RNA_NODUPS.CHIP_MULTIMARK.SORTED.subset.vcf'
with open(file_input) as vcf_handle:
    variant_set = vcf.read_vcf_variant_lines(vcf_handle, False)
    
for index in range(variant_set.length()):
        var_element = variant_set.seq[index]    # XXX: WARNING: changes made to element not saved

for index in range(variant_set.length()):
        var_element = variant_set.seq[index]        
        
        # 1. get reference sequence
        var_element.get_surround_seq(wing_l, fa_ind, force_ref_match)
        # 2. compute reverse complement
        var_element.assign_rev_complement()
        # 3. compute int version (faster to process as int)
        var_element.assign_int_versions()
        
        ref_seq = var_element.return_full_ref_seq_int(wing_l)
        var_seq = var_element.return_full_var_seq_int(wing_l)
        print("\tref int: " + format(ref_seq) +
              "\n\tvar int: " + format(var_seq))
        print("start motif_match_int")
        plusmatch = motif_set.motif_match_int(bp, ref_seq, var_seq, wing_l)
        
        print('## Positive Matches ##')
        for match in plusmatch:
            print( match.name + " vs=" + str(round(match.var_score, 4)) + 
                  " rs = " + str(round(match.ref_score, 4)) )
        
        
        # 6. Calculate motif matches to reverse complement
        ref_seq_rc = var_element.return_full_ref_seq_reverse_complement_int(wing_l)
        var_seq_rc = var_element.return_full_var_seq_reverse_complement_int(wing_l)
        print("\tref rc int: " + format(ref_seq_rc) +
              "\n\tvar rc int: " + format(var_seq_rc))
        print("start motif_match_int reverse complement")
        minusmatch = motif_set.motif_match_int(bp, ref_seq_rc, var_seq_rc, wing_l)
        
        print('## Reverse Complement Matches ##')
        for match in minusmatch:
            print( match.name + " vs=" + str(round(match.var_score, 4)) + 
                  " rs = " + str(round(match.ref_score, 4)) )<jupyter_output>pulling reference sequence for (chr1, 145039922, 1, 50, Fasta("../../data/genome_reference/reference_genome_hg19.fa"))
	145039872:145039973
	ref int: [3, 2, 3, 2, 3, 3, 2, 3, 2, 2, 2, 3, 3, 3, 3, 2, 4, 3, 4, 4, 4, 2, 3, 4, 3, 1, 3, 2, 4, 3, 3, 1, 3, 2, 2, 1, 3, 2, 3, 3, 2, 4, 3, 2, 2, 2, 3, 3, 3, 2, 3, 3, 2, 3, 3, 4, 3, 2, 3, 2, 3, 2, 1, 2, 4, 2, 4, 2, 2, 1, 3, 3, 2, 4, 3, 1, 3, 1, 2, 1, 2, 3, 1, 2, 4, 3, 3, 2, 4, 3, 3, 2, 1, 2, 3, 1, 3, 4, 4, 3, 2]
	var int: [3, 2, 3, 2, 3, 3, 2, 3, 2, 2, 2, 3, 3, 3, 3, 2, 4, 3, 4, 4, 4, 2, 3, 4, 3, 1, 3, 2, 4, 3, 3, 1, 3, 2, 2, 1, 3, 2, 3, 3, 2, 4, 3, 2, 2, 2, 3, 3, 3, 2, 2, 3, 2, 3, 3, 4, 3, 2, 3, 2, 3, 2, 1, 2, 4, 2, 4, 2, 2, 1, 3, 3, 2, 4, 3, 1, 3, 1, 2, 1, 2, 3, 1, 2, 4, 3, 3, 2, 4, 3, 3, 2, 1, 2, 3, 1, 3, 4, 4, 3, 2]
start motif_match_int
## Positive Matches ##
RFX5 vs=5.4636 rs = 8.4653
	ref rc int: [3, 2, 3, 2, 3, 3, 2, 3, 2, 2, 2, 3, 3, 3, 3, 2, 4, 3, 4, 4, 4, 2, 3, 4, 3, 1, 3, 2, 4, 3, 3, 1, 3, 2, 2, 1, 3, 2, 3, 3, 2, 4, 3, 2, 2, 2, 3, 3, 3,[...]
{'fraction_non_alphanumeric': 0.14490489754663236, 'fraction_numerical': 0.05494808416796839, 'mean_word_length': 4.275811924381968, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '428107', 'license': 'MIT', 'n_tokens_mistral': 4878, 'n_tokens_neox': 4257, 'n_words': 931, 'path': 'venusar/debugging_motifs.ipynb', 'repo_name': 'j-andrews7/VENUSAR', 'size': '10883'}