<jupyter_start><jupyter_text>Preprocess somatic mutations file<jupyter_code>import os
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import glob
import re
import numpy as np<jupyter_output><empty_output><jupyter_text>Set the path to save output files (PATH) and folder with data files (DATAPATH)<jupyter_code>PATH = "D:/Projects/cTaG2.0"
DATAPATH = "D:/Projects/data/cTaG2.0"<jupyter_output><empty_output><jupyter_text>Uploading the maf file to generate annovar input fileUpload the maf file<jupyter_code>os.chdir(DATAPATH + "/GDC_TCGA/SNV/995c0111-d90b-4140-bee7-3845436c3b42")
fname = "TCGA.BRCA.mutect.995c0111-d90b-4140-bee7-3845436c3b42.DR-10.0.somatic.maf"
data_snp=pd.read_csv(fname, sep="\t", comment="#", low_memory=False)
data_snp.shape
print("Total number of samples = {}".format(len(data_snp.Tumor_Sample_Barcode.unique())))
data_snp.head()
for col in data_snp.columns:
    print(col, data_snp.loc[0,col])<jupyter_output>Hugo_Symbol USP24
Entrez_Gene_Id 23358
Center WUGSC
NCBI_Build GRCh38
Chromosome chr1
Start_Position 55159655
End_Position 55159655
Strand +
Variant_Classification Missense_Mutation
Variant_Type SNP
Reference_Allele T
Tumor_Seq_Allele1 T
Tumor_Seq_Allele2 C
dbSNP_RS rs150880897
dbSNP_Val_Status by1000G;byCluster;byFrequency
Tumor_Sample_Barcode TCGA-D8-A1XY-01A-11D-A14K-09
Matched_Norm_Sample_Barcode TCGA-D8-A1XY-10A-01D-A14K-09
Match_Norm_Seq_Allele1 nan
Match_Norm_Seq_Allele2 nan
Tumor_Validation_Allele1 nan
Tumor_Validation_Allele2 nan
Match_Norm_Validation_Allele1 nan
Match_Norm_Validation_Allele2 nan
Verification_Status nan
Validation_Status nan
Mutation_Status Somatic
Sequencing_Phase nan
Sequence_Source nan
Validation_Method nan
Score nan
BAM_File nan
Sequencer Illumina HiSeq 2000
Tumor_Sample_UUID edb6d161-8f50-4c11-8246-487c4ea9a55d
Matched_Norm_Sample_UUID 8dea96d9-5017-4872-a84e-33bfd2f37b7a
HGVSc c.1024A>G
HGVSp p.Ile342Val
HGVSp_Short p.I342V
Transcript_ID ENST00000294383
[...]<jupyter_text>Generate input file for ANNOVAR<jupyter_code>data_annovar = data_snp.loc[:,["Chromosome", "Start_Position", "End_Position", "Reference_Allele", "Tumor_Seq_Allele2"]]
data_annovar.Chromosome = [x[3:] for x in data_annovar.Chromosome]
data_annovar.head()
# check end positions
temp=[]
for s,e,r,a in zip(data_annovar.Start_Position, data_annovar.End_Position, data_annovar.Reference_Allele, data_annovar.Tumor_Seq_Allele2):
    if s == e and len(r) == 1:
        temp.append(True)
    elif e-s == len(r) - 1:
        temp.append(True)
    else:
        temp.append(False)
temp.count(False)
        

len(data_annovar[data_annovar.Tumor_Seq_Allele2 == "-"])
os.chdir(DATAPATH + "/GDC_TCGA/SNV")
fname="brca_snp.avinput"
data_annovar.to_csv(fname, sep="\t", header=False, index=False)<jupyter_output><empty_output><jupyter_text>Run ANNOVARrun shell command to run annovar using the generated file as inputGo to folder containing annovar:cd /data/Priyanka/softwares/annovar/command:table_annovar.pl /data/malvika/temp/brca_snp.avinput /data/malvika/humandb/ -buildver hg38 -out brca_anno.hg38_multianno.csv -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f -nastring . -csvout -polish -xref example/gene_xref.txtCopy the resultant file to relevant folder  Concatenate the output file from ANNOVAR and maf fileLoad ANNOVAR output file<jupyter_code>os.chdir(DATAPATH + "/GDC_TCGA/SNV")
fname = "brca_anno.hg38_multiannov.hg38_multianno.csv"
data_anno=pd.read_csv(fname, sep=",", comment="#").drop(0).replace(".", np.nan).reset_index()
data_all = pd.concat([data_snp, data_anno], axis=1)<jupyter_output><empty_output><jupyter_text>Discard hypermutated samples<jupyter_code>max_mut = 2000
mut_count = data_all.groupby(["Tumor_Sample_Barcode"]).Tumor_Sample_Barcode.count()
data_all = data_all[~data_all.Tumor_Sample_Barcode.isin(mut_count[mut_count>2000].index)]
print("Total number of samples (after filtering) = {}".format(len(data_all.Tumor_Sample_Barcode.unique())))<jupyter_output>Total number of samples (after filtering) = 984<jupyter_text>Covert the old gene symbols into relevant genesLoad hugo to ensemble id file and check for gene symbols not recognised.Take unique list of genes and save to file.<jupyter_code># Load hugo to ensemble id mapping data
os.chdir("D:/Projects/data/cTaG2.0")
fname = "hugo2ensmbl.txt"
data_map = pd.read_csv(fname, sep="\t", header=0, index_col=1)

os.chdir(DATAPATH + "/GDC_TCGA/SNV")
fname = "convertgenes.txt"
pd.Series(list(set([x for x in data_all.Hugo_Symbol if x not in data_map.index]))).to_csv(fname, index=False, header=False)<jupyter_output><empty_output><jupyter_text>Use HGNC converterhttps://www.genenames.org/tools/multi-symbol-checker/If the converter throws error or takes entire gene list as one entry, delete and re-enter first new line in the uploaded genelist and try again.Download the output  as csv and open with notepad and delete first line from output csv file (sep=,).The output file is now ready for use.<jupyter_code>os.chdir(DATAPATH + "/GDC_TCGA/SNV")
fname = "hgnc-symbol-check.csv"
gene_map = pd.read_csv(fname, index_col=0, header=0)
gene_map = gene_map[gene_map["Match type"] != "Unmatched"]
gene_map.head()
# Update the gene symbols in dataframe
for idx, gene_sym, gene_ensg in zip(data_all.index, data_all["Hugo_Symbol"], data_all["Gene"]):
    if gene_sym in gene_map.index:
        if type(gene_map.loc[gene_sym, "Approved symbol"]) is not str:
            data_all.loc[idx, "Hugo_Symbol"] = (data_map[data_map["Ensembl gene ID"] == gene_ensg].index[0])
        else:
            data_all.loc[idx, "Hugo_Symbol"] = gene_map.loc[gene_sym, "Approved symbol"]
        
len(set([x for x in data_all.Hugo_Symbol if x not in data_map.index]))<jupyter_output><empty_output><jupyter_text>Keep relevant columns<jupyter_code>cols = ['Hugo_Symbol', "Gene", "Tumor_Sample_Barcode", "Chr", "Start", "End", "Ref",
        "Alt",'Variant_Classification', 'DOMAINS',
        'SIFT_score', 'SIFT_pred', 'Polyphen2_HDIV_score',
        'Polyphen2_HDIV_pred', 'Polyphen2_HVAR_score', 'Polyphen2_HVAR_pred',
        'LRT_score', 'LRT_pred', 'MutationTaster_score', 'MutationTaster_pred',
        'MutationAssessor_score', 'MutationAssessor_pred', 'FATHMM_score',
        'FATHMM_pred', 'PROVEAN_score', 'PROVEAN_pred', 'VEST3_score',
        'CADD_raw', 'CADD_phred', 'DANN_score', 'fathmm-MKL_coding_score',
        'fathmm-MKL_coding_pred', 'MetaSVM_score', 'MetaSVM_pred',
        'MetaLR_score', 'MetaLR_pred', 'integrated_fitCons_score',
        'integrated_confidence_value', 'GERP++_RS',
        'phyloP7way_vertebrate', 'phyloP20way_mammalian',
        'phastCons7way_vertebrate', 'phastCons20way_mammalian',
        'SiPhy_29way_logOdds']
data = data_all[cols]<jupyter_output><empty_output><jupyter_text>Save file<jupyter_code>os.chdir(DATAPATH + "/GDC_TCGA/SNV")
fname="brca_snv.tsv"
data.to_csv(fname, sep="\t", header=True, index=False)
data.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11547310900201323, 'fraction_numerical': 0.03508771929824561, 'mean_word_length': 6.336497890295359, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 41, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '71567', 'license': 'CC0-1.0', 'n_tokens_mistral': 2908, 'n_tokens_neox': 2700, 'n_words': 504, 'path': 'code/preprocess_SNV.ipynb', 'repo_name': 'RamanLab/PIVOT', 'size': '6954'}