from google.colab import auth
auth.authenticate_user()
print('Authorized')

"""
Explanation: Bam slicing in a cloud hosted jupyter notebook.
Welcome to the ‘Query of the Month’. This is part of our collection of new and interesting queries to demonstrate the powerful combination of BigData from the NCI cancer programs like TCGA, and BigQuery from Google.
Please let us know if you have an idea or a suggestion for our next QotM!
Query of the Month is produced by:
David L Gibbs (david.gibbs ( ~ at ~ ) systemsbiology ( ~ dot ~ ) org)
Kawther Abdilleh (kawther.abdilleh ( ~ at ~ ) gdit (~ dot ~) com)
Sheila M Reynolds (sheila.reynolds ( ~ at ~ ) systemsbiology ( ~ dot ~ ) org)
In this notebook, using the Pysam package, we demonstrate how to slice bam files stored in GCS buckets.
Here, you'll learn how to: 

How to invoke bash commands within a Jupyter environment. 
How to install packages/programs within a Jupyter environment
How to use available BigQuery tables within ISB-CGC to query and identify Google Cloud Storage bucket locations for BAM files of interest
How to use PySam to slice BAM files
How to save slices in your bucket and retrieve them
Brief example of working with reads

Let's authenticate ourselves
End of explanation
"""


import os
os.environ['HTSLIB_CONFIGURE_OPTIONS'] = "--enable-gcs"

"""
Explanation: First, let's prep to install Pysam and the HTSlib
Pysam is a python wrapper around samtools, and samtools uses the HTSlib (http://www.htslib.org/). So we need to make sure we have the necessary libraries to compile HTSlib and samtools. The compilation is needed to activate the ability to read from google cloud buckets.
End of explanation
"""



!ls -lha

!sudo apt-get install autoconf automake make gcc perl zlib1g-dev libbz2-dev liblzma-dev libcurl4-openssl-dev libssl-dev

!pip3 install pysam -v --force-reinstall --no-binary :all:

# Without forcing the compilation, we get error 
# [Errno 93] could not open alignment file '...': Protocol not supported

import pysam

"""
Explanation: We can invoke bash commands to see what was downloaded into our current working directory. Bash commands can be invoked by putting an exclamation point (!) before the command.
End of explanation
"""


# First, we need to set our project. Replace the assignment below
# with your project ID.
# project_id = 'isb-cgc-02-0001'
#!gcloud config set project {project_id}

#import os
#os.environ['GCS_OAUTH_TOKEN'] = "gcloud auth application-default print-access-token"


"""
Explanation: First, we need to set our project. Replace the assignment below with your project ID.
End of explanation
"""


%%bigquery --project isb-cgc-02-0001 df
SELECT * FROM `isb-cgc.TCGA_hg19_data_v0.tcga_metadata_data_hg19_18jul`
where 
data_format = 'BAM' 
AND disease_code = 'OV' 
AND experimental_strategy = "WGS" 
AND platform = 'ABI SOLiD'
LIMIT 5


"""
Explanation: Now that we have Pysam installed, let's write an SQL query to locate BAM files in Google Cloud Storage Buckets. 
In the query below, we are looking to identify the Google Cloud Storage bucket locations for TCGA Ovarian Cancer BAMs obtained via whole genome sequencing (WGS) generated using the SOLiD sequencing system
End of explanation
"""



samfile = pysam.AlignmentFile('gs://isb-ccle-open/gdc/0a109993-2d5b-4251-bcab-9da4a611f2b1/C836.Calu-3.2.bam', "rb")
for read in samfile.fetch('7', 140453130, 140453135):
  print(read)

samfile.close()

"""
Explanation: Now using the following Pysam command,  let's read a bam file from GCS and slice out a section of the bam using the fetch function. For the purposes of the BAM slicing exercise, we will use an open-access CCLE BAM File open-access BAM file.  CCLE open access BAM files are stored here
End of explanation
"""


samfile = pysam.AlignmentFile('gs://isb-ccle-open/gdc/0a109993-2d5b-4251-bcab-9da4a611f2b1/C836.Calu-3.2.bam', "rb")
fetchedreads = pysam.AlignmentFile("test.bam", "wb", template=samfile)
for read in samfile.fetch('7', 140453130, 140453135):
  fetchedreads.write(read)

fetchedreads.close()
samfile.close()

"""
Explanation: The output from the above command is a tab-delimited human readable table of a slice of the BAM file. This table gives us information on reads that mapped to the region that we "extracted" from chromosome 7 between the coordinates of 140453130 and 140453135.
Now, let's suppose you would like to save those reads to your own Google cloud storage bucket...
End of explanation
"""


!ls -lha

#if you don't already have a google cloud storage bucket, you can make one using the following command:
#The mb command creates a new bucket. 
#gsutil mb gs://your_bucket

#to see what's in the bucket..
#!gsutil ls gs://your_bucket/



# then we can copy over the file

!gsutil cp gs://bam_bucket_1/test.bam test_dl.bam

# and it made it?

!gsutil ls gs://bam_bucket_1/

"""
Explanation: Let's see if we saved it?
End of explanation
"""


newsamfile = pysam.AlignmentFile('gs://bam_bucket_1/test.bam', 'rb')
for r in newsamfile.fetch(until_eof=True):
  print(r)
#  
#  
# No. But maybe soon. #
  

"""
Explanation: Now, can we read it back?!?
End of explanation
"""


!gsutil ls gs://bam_bucket_1/

!gsutil cp gs://bam_bucket_1/test.bam test_dl.bam

"""
Explanation: Let's move our slice back to this instance.
End of explanation
"""


import numpy as np

# first we'll open our bam-slice
dlsamfile = pysam.AlignmentFile('test_dl.bam', 'rb')

# and we'll save the read quality scores in a list
quality = []
for read in dlsamfile:
  quality.append(read.mapping_quality)
  
# then we can compute statistics on them
print("Average quality score")
print(np.mean(quality))

# again open our bam-slice
dlsamfile = pysam.AlignmentFile('test_dl.bam', 'rb')

# here, we'll extract the sequences to process
seqs = []
for read in dlsamfile:
  seqs.append(read.query_sequence)
  

# let's count the number of times nucleotide bases are read.
baseCount = dict()
baseCount['A'] = 0
baseCount['C'] = 0
baseCount['G'] = 0
baseCount['T'] = 0
baseCount['N'] = 0
totalCount = 0

for si in seqs:
  for base in si:
    baseCount[base] += 1
    totalCount +=1

for bi in ['A', 'G', 'C', 'T']:
  baseCount[bi] /= totalCount
    
print("Percent bases")
print(baseCount)

"""
Explanation: Now we're ready to work with our bam-slice!
Very brief examples of working with reads.
The Alignment file is read as a pysam.AlignedSegment, which is a python class. The methods and class variables can be found here: https://pysam.readthedocs.io/en/latest/api.html#pysam.AlignedSegment
End of explanation
"""


