<jupyter_start><jupyter_text>From ChiP to Core Regulatory Circuitry.<NAME>@jkobject[www.jkobject.com](https://www.jkobject.com)in partnership with the Broad Institutethanks to young Lab @ MIT<jupyter_code>from Bio import SeqIO
import os
path_to_rose = "/Users/jeremie/Documents/logiciels/rose"
mypath= "datatosend"
os.chdir(path_to_rose)
import ROSE_utils
a = ["/Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC"] #! pwd
a = a[0]
%load_ext autoreload
%autoreload 2<jupyter_output><empty_output><jupyter_text>you might need to install an extention to curl to download from an sftp server. otherwise reinstall the latest version of curl with libsshon mac: `brew install curl --with-libssh2=/usr/local`after installing [conda]( ```conda config --add channels defaultsconda config --add channels biocondaconda config --add channels conda-forge```then`conda config install bowtie2 && pip install MACS2`then `git clone` and provide path to [ROSE](https://bitbucket.org/young_computation/rose/src) and [CRCMapper](https://bitbucket.org/young_computation/crcmapper)for CRC mapper you need [meme](http://meme-suite.org/doc/install.html?man_type=web)Not every version will work as they do not contain the same functions unfortunatly. I have found version **4.12.0** to work.<jupyter_code>mypath= "datatosend"
! curl -k "sftp://bioinfodata.wistar.upenn.edu/datatosend.tar.gz" --user [NAME] > $mypath.tar.gz<jupyter_output><empty_output><jupyter_text>the call might stop for many reasons. in this case just run this function below to continue it from where you started as the -R option is not working well with curl and sftp<jupyter_code>def notfinished():
    s = os.path.getsize(mypath+'.tar.gz')
    ! curl -C $s -k "sftp://bioinfodata.wistar.upenn.edu/datatosend.tar.gz" --user [NAME] > $mypath1.tar.gz
    ! cat $mypath.tar.gz $mypath1.tar.gz > $mypath2.tar.gz
    ! rm $mypath.tar.gz $mypath1.tar.gz
    ! mv $mypath2.tar.gz $mypath.tar.gz<jupyter_output><empty_output><jupyter_text>here we prepare the libs and the alignments with the hg19/37 version of the genome as it is the best one currently working with CRCmapper<jupyter_code>! gzip -d $mypath.tar.gz
! curl ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.dna.toplevel.fa.gz > reference.fa.gz
! gzip -d reference.fa.gz
! mkdir reference
os.chdir("reference")
! wget --timestamping 'ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y}.fa.gz'
! gzip -d *.fa.gz<jupyter_output><empty_output><jupyter_text>some visualization<jupyter_code>for i, rec in enumerate(SeqIO.parse(mypath+"/54-mp457-PDX34-H3K27ac-r2_S7_R1.fastq", "fastq")):
    print(rec)
    if i > 10:
        break<jupyter_output><empty_output><jupyter_text>here we can see the -logp that the base call was incorrect (called phred quality/sanger quality...)<jupyter_code>for i, rec in enumerate(SeqIO.parse(mypath+"/54-mp457-PDX34-H3K27ac-r2_S7_R1.fastq", "fastq")):
    print(rec.seq)
    if i > 10:
        break
! MACS2 callpeak -h
! bowtie2 -h<jupyter_output><empty_output><jupyter_text>Let's build the indexes with bowtie2<jupyter_code>! bowtie2-build $mypath/reference.fa index --threads 8
#maybe try to output bampe files next time (and with higher sensibility) for macs2 to read the paired end information<jupyter_output><empty_output><jupyter_text>we then produce the alignmentshere we use the paired end alignments<jupyter_code>! bowtie2 -x index -1 $mypath/54-mp457-PDX34-H3K27ac-r2_S7_R1.fastq.gz -2 $mypath/54-mp457-PDX34-H3K27ac-r2_S7_R2.fastq.gz -S $mypath/mapped.sam --threads 8 -t<jupyter_output><empty_output><jupyter_text>we can then use macs2 for the peak calling<jupyter_code>! MACS2 callpeak -f BAMPE -t $mypath/mapped.sam --outdir "peak_called_bampe"
#<jupyter_output><empty_output><jupyter_text>Some visualization...<jupyter_code>! Rscript peak_called/NA_model.r
from IPython.display import IFrame
IFrame("./peak_called/NA_model.pdf", width=600, height=300)<jupyter_output><empty_output><jupyter_text>finally we need some file transformation for rose to be able to process everything wellwe transform the outputed bed files into gff for rose to be able to process<jupyter_code>import ROSE_utils
ROSE_utils.bedToGFF(bed="/Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC/peak_called/NA_summits.bed",output="/Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC/peak_called/NA_summits.gff")<jupyter_output><empty_output><jupyter_text>we then sort and index the mapped reads<jupyter_code>! samtools sort --threads 6 $mypath/mapped.sam > $mypath/mapped_sorted.bam<jupyter_output><empty_output><jupyter_text>we erase the now useless file<jupyter_code>! rm $mypath/mapped.sam<jupyter_output><empty_output><jupyter_text>we sort the reads<jupyter_code>! samtools index -@ 7 $mypath/mapped_sorted.bam<jupyter_output><empty_output><jupyter_text>finally we have to switch folders and then call roserose will now stitch the enhancers to find the superenhancers in this cell line<jupyter_code>import os
a = ['/Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC'] #! pwd
a = a[0]
path_to_rose = "/Users/jeremie/Documents/logiciels/rose"
os.chdir(path_to_rose)
! python ROSE_main.py -g hg19 -i $a/peak_called/NA_summits.gff -r $a/$mypath/mapped_sorted.bam -o $a/superenhancers<jupyter_output><empty_output><jupyter_text>We can now create a map of the core regulatory circuitry of this cell line! they were two lines to change in crcmapper to make it work...```pythondef createExpressionDict(annotationFile, projectFolder, projectName, refseqToNameDict, expressionTable):```-----line 98----and ```pythonfor line in expressionTable[1:]:    trid = line[0]    geneName = refseqToNameDict[trid]    if len(expressionTable[1]) == 3:   when expressionTable is an output from bamToGFF.py        exp = float(line[2])    else:   when expressionTable is passed as an option (2 columns)        exp = float(line[1])```----line 120 creating syntax error----and```pythonself._chr = self.__chrDict[chr] if self.__chrDict[chr][:3] == 'chr' else 'chr' + self.__chrDict[chr]```----in utils.py line 44to accomodate the potential variation in name in bam output files (chr or ) ---and```python382:done = False    for line in fimoTable[1:]:        source = motifDatabaseDict[line[0]]        if not done:            for p, s in enumerate(line):                if '|' in s:                    pos = p                    done = True                    break        region = line[pos].split('|')403:if (region[1], int(region[2]) + int(line[pos + 1]), int(region[2]) + int(line[pos + 2])) not in motifDictSE[source]:405:motifDictSE[source].append((region[1], int(region[2]) + int(line[pos + 1]), int(region[2]) + int(line[pos + 1])))```---to accomodate for the potential different positions of the fimo table output according to different `meme` version---<jupyter_code>path_to_crcmapper = "/Users/jeremie/Documents/logiciels/crcmapper"
os.chdir(path_to_crcmapper)
from CRCmapper import findCRC
print a, mypath
findCRC(
    fasta=a+"/"+mypath+"/reference/", 
    genome="HG19", 
    enhancers=a+"/superenhancers/NA_summits_AllEnhancers.table.txt", 
    bam=a+"/"+mypath+"/mapped_sorted.bam", 
    subpeaks=a+"/peak_called/NA_summits.bed", 
    name="CRC", 
    output=a+"/crcmapper", 
    activity=a+"/crcmappermatrix.gff",
    expCutoff=33, extension=500, Enumber="supers")
! mkdir CRC
! mv crcmapper* CRC/<jupyter_output>/Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC datatosend
IDENTIFY EXPRESSED GENES
python ./bamToGFF.py -r  -d  -o /Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC/crcmappermatrix.gff -m 1 -f 0 -e 200  -i /Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC/crcmapperCRC_TSS.gff -b /Users/jeremie/Documents/Projects/BroadInstitute/CHiP2CRC/datatosend/mapped_sorted.bam
CREATING SUPER-ENHANCER LOCUS COLLECTION
CREATING EXPRESSION DICTIONARY
Expression cutoff: 0.0
FINDING CANDIDATE TFs
Number of candidate TFs: 27
CREATE CANDIDATE TFs AND SE TABLE
MAKE FASTA
MOTIF SEARCH
Number of annotated candidate TFs that have motifs: 14
fimo --motif 'JASPAR2014.MA0102.3' --motif 'Transfac.V$CFOS_Q6' --motif 'JASPAR2014.MA0476.1' --motif 'Transfac.V$JUNB_Q6' --motif 'JASPAR2014.MA0490.1' --motif 'Transfac.V$KLF3_Q3' --motif 'Transfac.V$KLF7_04' --motif 'Transfac.V$KLF7_03' --motif 'Uniprobe.UP00093_1' --motif 'Uniprobe.UP00093_2' --motif 'Transfac.V$MEF2D_01' --motif 'Jolma2013.M[...]<jupyter_text>an automated version<jupyter_code>mypath = "datatosend"
refs = [30,31,32,33,34]
files = {}
path_to_rose = "/Users/jeremie/Documents/logiciels/rose"
path_to_crcmapper = "/Users/jeremie/Documents/logiciels/crcmapper"
a = ! pwd
a = a[0]
os.chdir(path_to_rose)
import ROSE_utils
os.chdir(a)
for ref in refs
    for r in os.listdir(mypath):
        if r[:2] in refs:
            lane = "R1" if "_R1" in r else "R2"
            files[lane] = r
            ! bowtie2 -x index -1 $mypath/$files["R1"] -2 $mypath/$files["R2"] -S $mypath/mapped.sam --threads 8 -t
            ! MACS2 callpeak -t $mypath/mapped.sam --outdir peak_called_$ref
            ! Rscript peak_called_$ref/NA_model.r
            ! samtools sort --threads 6 $mypath/mapped.sam > $mypath/mapped_sorted_$ref.bam 
            ! rm $mypath/mapped.sam
            ! samtools index -@ 7 $mypath/mapped_sorted_$ref.bam 
            os.chdir(path_to_rose)
            ROSE_utils.bedToGFF(bed=a + "/peak_called_"+ref+"/NA_summits.bed",output=a + "/peak_called_"+ ref +"/NA_summits.gff")
            ! python ROSE_main.py -g hg19 -i $a/peak_called_$ref/NA_summits.gff -r $a/$mypath/mapped_sorted_$ref.bam -o $a/superenhancers_$ref
            os.chdir(path_to_crcmapper)
            ! python CRCmapper.py -e $a/superenhancers_$ref/NA_summits_AllEnhancers.table.txt -b $a/$mypath/mapped_sorted_$ref.bam -g hg19 -f $a/$mypath/reference  -s $a/peak_called_$ref/NA_summits.bed -n CRC -o $a/crcmapper_$ref
            os.chdir(a)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12781954887218044, 'fraction_numerical': 0.028471177944862154, 'mean_word_length': 4.545302946081156, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 77, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '481314', 'license': 'Apache-2.0', 'n_tokens_mistral': 4005, 'n_tokens_neox': 3718, 'n_words': 802, 'path': 'from ChiP to Core Regulatory Circuitry.ipynb', 'repo_name': 'jkobject/from-ChiP-to-Core-Regulatory-Circuitry1', 'size': '9983'}