<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 sequana_coverage test case example (fungus)
This notebook creates the BED file S_pombe.filtered.bed provided in 
- https
Step1:  Download FastQ files (1.6Gb)
Step2:  Download reference and annotation files
Step3:  The Reference must be altered to rename the header so that they 
agree with the genbank
Step4:  Sequana_coverage analysis (using the library)
Step5:  Chromosome selection
This data set contains several chromosomes. We can select them one by one and anlyse them as follows
Step6:  Sequana_coverage analysis (using the standalone)
<ASSISTANT_TASK:>
Python Code:
%pylab inline
matplotlib.rcParams['figure.figsize'] = [10,7]
Explanation: sequana_coverage test case example (fungus)
This notebook creates the BED file S_pombe.filtered.bed provided in 
- https://github.com/sequana/resources/tree/master/coverage and
- https://www.synapse.org/#!Synapse:syn10638358/wiki/465309
genome length: 5.5Mb
It also shows the ability of sequana_coverage tool to handle multi-chromosome input data sets
WARNING: To create the input BED file, you will need an account on synapse to get the FastQ files.
You can skip the steps that build the BED file and download it directly from github:
wget //github.com/sequana/resources/raw/master/coverage/S_pombe.filtered.bed.bz2
bunzip2 S_Pombe.filtered.bed
and jump to the Sequana_coverage analysis (using the library) section directly.
Otherwise, we first download 2 FastQ files from Synapse, its reference genome and its genbank annotation. Then, we use   BWA to map reads into a BAM file. The BAM file itself is converted to a BED, which is going to be one input file to our analysis. Finally, we use the coverage tool from Sequana project (i) with the standalone (sequana_coverage) and (ii) the Python library to analyse the BED file.
Versions used:
- sequana 0.7.0
- bwa mem 0.7.15
- bedtools 2.26.0
- samtools 1.5
- synapseclient 1.7.2
End of explanation
import synapseclient
l = synapseclient.login()
l.get("syn10641621", downloadLocation=".", ifcollision="overwrite.local")
l.get("syn10641896", downloadLocation=".", ifcollision="overwrite.local")
Explanation: Download FastQ files (1.6Gb)
End of explanation
!sequana_coverage --download-reference CU329670 
!sequana_coverage --download-reference CU329671
!sequana_coverage --download-reference CU329672 
!sequana_coverage --download-reference X54421 
!sequana_coverage --download-genbank CU329670
!sequana_coverage --download-genbank CU329671
!sequana_coverage --download-genbank CU329672
!sequana_coverage --download-genbank X54421
!cat CU*gbk  X*gbk > S_pombe.gbk
!cat CU*.fa X*.fa> S_pombe.fa
Explanation: Download reference and annotation files
End of explanation
files = ['CU329670.fa', "CU329671.fa", "CU329672.fa", "X54421.fa"]
with open("S_pombe.fa", "w") as fout:
    for filename in files:
        with open(filename, "r") as fin  :
            for line in fin.readlines():
                if line.startswith(">"):
                    start, end = line.split(None, 1)
                    accession = start[1:].rsplit("|", 1)[1]
                    line = ">" + accession + " " + end
                fout.write(line)
# The mapping to obtain the sorted BAM file (uses BWA behing the scene)
!time sequana_mapping \
    --file1 M14-19_J29_01_TAAGGCGA-TATCCTCT_L002_R1_001.fastq.gz \
    --file2 M14-19_J29_01_TAAGGCGA-TATCCTCT_L002_R2_001.fastq.gz \
    --reference S_pombe.fa --thread 4
# Build the BED file (unfiltered)
! time samtools depth -d 30000 S_pombe.fa.sorted.bam -aa > S_pombe.bed
Explanation: The Reference must be altered to rename the header so that they 
agree with the genbank
End of explanation
%%time
from sequana import GenomeCov
# If chromosome length is >5Mb, we split the data. Here we now it is 5.5Mb, so let us
#slightly increase the chunksize. 
b = GenomeCov("S_pombe.bed", "S_pombe.gbk", chunksize=6000000, low_threshold=-4, high_threshold=4)
b.compute_gc_content("S_pombe.fa")
Explanation: Sequana_coverage analysis (using the library)
End of explanation
chrom = b.chr_list[0]
chrom.run(20001, circular=True)
chrom.plot_coverage()
_ = ylim([0, 400])
chrom.plot_rois(0,1000000)
chrom.plot_gc_vs_coverage(bins=[80, 60], Nlevels=6)
Explanation: Chromosome selection
This data set contains several chromosomes. We can select them one by one and anlyse them as follows:
End of explanation
!sequana_coverage --input S_pombe.filtered.bed --genbank S_pombe.gbk --reference S_pombe.fa
Explanation: Sequana_coverage analysis (using the standalone)
End of explanation
<END_TASK>