<jupyter_start><jupyter_text>get padded cds and transcript sequences for transcript set chosen with ORFquant<jupyter_code>import pandas as pd
import numpy as np
import math

def convertFastaToDict(fastaFile):
    '''
    converts a fasta file to a dict of {sequenceName:sequence}
    can take extra files in * args
    '''
    if isinstance(fastaFile, list):
        files = fastaFile
    else:
        files = [fastaFile]
    currentName = None
    currentSequence = None
    seqDict = {}
    for currentFile in files:
        if currentFile.endswith('.gz'):
            import gzip
            f = gzip.open(currentFile)
        else:
            f = open(currentFile)
        for line in f:
            if not line.strip() == '' and not line.startswith('#'):  # ignore empty lines and commented out lines
                if line.startswith('>'):  # > marks the start of a new sequence
                    if not currentName == None:  # after we've reached the firtst > line, we know what the sequence corresponds to
                        seqDict[currentName] = currentSequence
                    currentName = line.strip()[1:].split()[
                        0]  # i've noticed the gencode names have extraneous numbering after some whitespace. This doens't match the GTF files, so I'm removing it.
                    currentSequence = ''
                else:
                    currentSequence += line.strip()
        f.close()
    seqDict[currentName] = currentSequence
    return seqDict

def get_CDS_seq(transcript, pad_5p, pad_3p):
    transcript_seq = transcript_seqs[transcript].upper()
    info_parts = transcript.split('|')
    for info_part in info_parts:
        if info_part.startswith('CDS:'):
            #print info_part, transcript
            CDS_positions = [int(pos) for pos in info_part.split(':')[1].split('-')]
            CDS_seq = transcript_seq[max((CDS_positions[0]-1)-pad_5p, 0):CDS_positions[1]+pad_3p]
            return CDS_seq
#download transcript sequences from gencode
! curl -L ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/gencode.v35.pc_transcripts.fa.gz -o ../annotations/gencode.v35.pc_transcripts.fa.gz
! gunzip -c ../annotations/gencode.v35.pc_transcripts.fa.gz | cat - ../annotations/reporter/UGAC_reporter_transcript.fa > ../annotations/gencode.v35.pBZ105.pc_transcripts.fa
transcript_seqs = convertFastaToDict('../annotations/gencode.v35.pBZ105.pc_transcripts.fa')
def get_tx_from_file(file_name):
    txs = set()
    f = open(file_name)
    for line in f:
        txs.add(line.strip())
    f.close()
    return txs
def write_padded_CDS_seqs(transcript_file, outfile_name, pad_5p, pad_3p):
    transcript_subset = get_tx_from_file(transcript_file)
    outfile = open(outfile_name, 'w')
    for transcript in transcript_seqs.keys():
        if transcript.split('|')[0] in transcript_subset:
            cds_seq = get_CDS_seq(transcript, pad_5p, pad_3p)
            outfile.write('>%s\n%s\n' % (transcript, cds_seq))
    outfile.close()

def write_subset_tx_seqs(transcript_file, outfile_name):
    transcript_subset = get_tx_from_file(transcript_file)
    outfile = open(outfile_name, 'w')
    for transcript in transcript_seqs.keys():
        if transcript.split('|')[0] in transcript_subset:
            outfile.write('>%s\n%s\n' % (transcript, transcript_seqs[transcript]))
    outfile.close()
tx_file = '../1_orfquant/orfquant_tx_stops_collapsed.txt'
import os
write_padded_CDS_seqs(tx_file, os.path.join('../annotations/', 'orfquant_CDS_%d_%d.fa' % (15, 15)),15, 15)
write_subset_tx_seqs(tx_file, '../annotations/orfquant_tx.fa')<jupyter_output><empty_output><jupyter_text>generate salmon indexes<jupyter_code>#generate decoys list, which just contains a list of human chromosome names
! grep "^>" < ../annotations/genome/GRCh38.primary_assembly.genome.pBZ105.fa | cut -d " " -f 1 > decoys.txt
! sed -i.bak -e 's/>//g' decoys.txt
#now generate "gentrome" files, which concatenate the genome and transcriptome
!mkdir gentromes
import os
import subprocess
for tx_file in ['orfquant_CDS_15_15.fa', 'orfquant_tx.fa']:
    tx_path = os.path.join('../annotations/', tx_file)
    prefix = tx_file[:-3]
    gentrome_path = os.path.join('gentromes/', prefix+'.gentrome.fa.gz')
    cmd="cat {tx_path} ../annotations/genome/GRCh38.primary_assembly.genome.pBZ105.fa | gzip > {gentrome_path}".format(**locals())
    print(cmd)
    #subprocess.Popen(cmd, shell=True).wait()
    ! {cmd}
#now generate a salmon index from each gentrome
!mkdir indices
import os
import subprocess
tx_files = ['orfquant_CDS_15_15.fa', 'orfquant_tx.fa']
for tx_file in sorted(tx_files):
    for k in [13, 31]:
        prefix = tx_file[:-3]
        gentrome_path = os.path.join('gentromes/', prefix+'.gentrome.fa.gz')
        index_path = os.path.join('indices/', prefix)
        #use K=31 for RNAseq (~75bp reads, and 13 for profiling)
        cmd="salmon index -t {gentrome_path} -d decoys.txt t -p 40 -i {index_path}_k{k} -k {k} --gencode ".format(**locals())
        #pasted into terminal
        print(cmd)
        ! {cmd}
    #subprocess.Popen(cmd, shell=True).wait()<jupyter_output>mkdir: cannot create directory ‘indices’: File exists
salmon index -t gentromes/orfquant_CDS_15_15.gentrome.fa.gz -d decoys.txt t -p 40 -i indices/orfquant_CDS_15_15_k13 -k 13 --gencode 
Version Info: This is the most recent version of salmon.
[2020-11-24 21:39:01.728] [jLog] [info] building index
out : indices/orfquant_CDS_15_15_k13
[00m[2020-11-24 21:39:01.728] [puff::index::jointLog] [info] Running fixFasta
[00m
[Step 1 of 4] : counting k-mers

[35m[2020-11-24 21:40:18.465] [puff::index::jointLog] [warning] Removed 90 transcripts that were sequence duplicates of indexed transcripts.
[00m[35m[2020-11-24 21:40:18.465] [puff::index::jointLog] [warning] If you wish to retain duplicate transcripts, please use the `--keepDuplicates` flag
[00m[00m[2020-11-24 21:40:18.466] [puff::index::jointLog] [info] Replaced 151,122,963 non-ATCG nucleotides
[00m[00m[2020-11-24 21:40:18.466] [puff::index::jointLog] [info] Clipped poly-A tails from 3 transcripts
[00mwrote 24440 cleaned reference[...]
{'fraction_non_alphanumeric': 0.12681336593317033, 'fraction_numerical': 0.038793806030969845, 'mean_word_length': 3.469045884923525, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '175279', 'license': 'MIT', 'n_tokens_mistral': 2314, 'n_tokens_neox': 2053, 'n_words': 492, 'path': 'F4_profiling/2_make_salmon_indices/generate_indexes.ipynb', 'repo_name': 'borisz264/NMD_screen_2020', 'size': '6135'}