%pylab inline
matplotlib.rcParams['figure.figsize'] = [10,7]

"""
Explanation: sequana_coverage test case example (Bacteria Staphylococcus aureus)
This notebook creates the BED file provided in 
- https://github.com/sequana/resources/tree/master/coverage and
- https://www.synapse.org/#!Synapse:syn10638358/wiki/465309
First, we download 2 FastQ from EBI FTP, its reference genome and its genbank annotation. Then, we map the read onto the reference using BWA and obtain a BAM file. The BAM file is converted to a BED, which is going to be one input file for our analysis. Finally, we use the coverage tool from Sequana project (i) with the standalone (sequana_coverage).
At the bottom of the page, we also added an example of an application where we show that some SNPs could be discarded or assigned a negative score when they appear next or within a ROIs.
Versions used:
- sequana 0.7.0
- bwa mem 0.7.15
- bedtools 2.26.0
- samtools 1.5
End of explanation
"""


from bioservices import ENA
ena = ENA()
data = ena.get_data("FN433596", frmt="fasta")
with open("FN433596.fasta", "wb") as fh:
    fh.write(data)

"""
Explanation: Download the genbank and genome reference
We will use bioservices (ENA service) or you can download the file
from our github: https://github.com/sequana/resources/tree/master/coverage
End of explanation
"""


!wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR036/ERR036019/ERR036019_1.fastq.gz
!wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR036/ERR036019/ERR036019_2.fastq.gz

"""
Explanation: Download the fastq files
End of explanation
"""


!sequana_mapping --file1 ERR036019_1.fastq.gz --file2 ERR036019_2.fastq.gz --reference FN433596.fasta --thread 4

"""
Explanation: Map the reads
depending on your system this step may take a while (10-30mins)
End of explanation
"""


!bedtools genomecov -d -ibam FN433596.fasta.sorted.bam > FN433596.bed

"""
Explanation: Convert BAM to BED
End of explanation
"""


!samtools view -q 35 -o FN433596.filtered.bam FN433596.fasta.sorted.bam
!samtools depth -d 20000 FN433596.fasta.sorted.bam FN433596.filtered.bam -aa > FN433596.filtered.bed

"""
Explanation: Filter the BAM file to get a 4-columns BED file (optional)

first column, as usual contains chromosome name
second column contains the position
third column contains the coverage
fourth column contains the coverage for mapped reads with high quality
End of explanation
"""


from sequana.snpeff import download_fasta_and_genbank
download_fasta_and_genbank("FN433596", "FN433596", fasta=False)

"""
Explanation: get a genbank (for the standalone sequana_coverage)
End of explanation
"""


from sequana import GenomeCov
# load the BED and genbank file
b = GenomeCov("FN433596.bed", "FN433596.gbk")
b.compute_gc_content("FN433596.fasta")

# load the first chromosome and compute running median and get z-scores
chromosome = b.chr_list[0]
chromosome.run(30001, circular=True)

# 5Mb data points. So, this may take a while. Note that only 1Mb are used for plotting.
chromosome.plot_coverage()

"""
Explanation: Sequana coverage (library)
End of explanation
"""


#!sequana_coverage --input FN433596.bed --reference FN433596.fasta --genbank FN433596.gbk -w 80001
!sequana_coverage --input FN433596.bed --reference FN433596.fasta  -w 80001 -o

"""
Explanation: You can also use the sequana_coverage standalone application
End of explanation
"""


!cp report/coverage_reports/*/rois.csv .

"""
Explanation: You can now open the report ./report/multiqc_report.html
End of explanation
"""


# The file ERR036019.filter.vcf was generated with sequana pipeline variantcalling
# version 0.6.2 using sequanix interface. The config file and pipeline 
# can be found in ./pipeline directory

# The rois.csv file was generated by the sequana_coverage command here above
# in 
import pandas as pd
roi = pd.read_csv("rois.csv")

from sequana.freebayes_vcf_filter import VCF_freebayes as VCF
vcf = VCF("ERR036019.filter.vcf")
variants = [v for v in vcf]
positions = [v.POS for v in variants]
quality = [v.QUAL for v in variants]  # freebayes score
M = max(quality)
scores = [x/M*100 for x in quality] # normalise
depths = [chromosome.df['cov'].iloc[pos] for pos in positions]

def plot_variant(x1=526000, x2=538000):
    clf()
    chromosome.plot_coverage()
    low = roi[roi.mean_cov<roi.mean_rm]
    high = roi[roi.mean_cov>roi.mean_rm]

    chr = chromosome
    highT = (chr.thresholds.high * chr.best_gaussian["sigma"] +
         chr.best_gaussian["mu"]) * chr.df["rm"]
    lowT = (chr.thresholds.low * chr.best_gaussian["sigma"] +
        chr.best_gaussian["mu"]) * chr.df["rm"]

    for k,v in low.iterrows():        
        Y1 = chr.df['cov'].iloc[v.start:v.end]
        Y2 = lowT.iloc[v.start:v.end] 
        Y1 = Y1.combine(Y2, max) *0
        if v.start > x1 and v.end<x2:
            fill_between(range(v.start, v.end), Y1, Y2, alpha=0.6, color="blue")
    
    for k,v in high.iterrows():
        Y1 = highT.iloc[v.start:v.end]
        Y2 = chr.df['cov'].iloc[v.start:v.end]
        Y2 = Y2.combine(Y1,max)
        if v.start > x1 and v.end<x2:
            fill_between(range(v.start, v.end), y1=Y1,y2=Y2 ,alpha=0.6, color="orange")

    scatter(positions, depths, c=scores, s=150)
    colorbar()
    xlim([x1,x2])

plot_variant()

"""
Explanation: Application: SNP annotation
In the directory, we provide a file called ERR036019.filter.vcf that contains 
a VCF file generated with the Variant Calling pipeline described here
http://sequana.readthedocs.io/en/master/pipeline_variant_calling.html
Here below, we read the VCF, the coverage and the ROIs found by the sequana_coverage tool. For a few examples, we plot the coverage, the ROIS (blue areas) and the SNPs found. We show that SNPs we low quality also coincide with ROIs. They may be further discarded and confirmed to be poor quality SNPs.
End of explanation
"""


plot_variant(1913000,1923000)

plot_variant(947000, 952000)

"""
Explanation: Here, the color in the RHS indicates the quality of the SNPs
End of explanation
"""


