<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Notebook Analyzes Duplication in Motif Defintions
Step1:  Import Motif Definitions
Step2:  Convert to DataFrame for Analysis of Duplication
Step3:  Move On
<ASSISTANT_TASK:>
Python 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)  
Explanation: Notebook Analyzes Duplication in Motif Defintions
End of explanation
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))
Explanation: Import Motif Definitions
End of explanation
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']
Explanation: Convert to DataFrame for Analysis of Duplication
End of explanation
# 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)) )
        
        
Explanation: 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)
End of explanation
<END_TASK>