Merge pull request #20 from luslab/dev

## Pipeline summary

1. Download FastQ files via SRA, ENA or GEO ids and auto-create input samplesheet ([`ENA FTP`](https://ena-docs.readthedocs.io/en/latest/retrieval/file-download.html); *if required*)

2. Merge re-sequenced FastQ files ([`cat`](http://www.linfo.org/cat.html))

3. Read QC ([`FastQC`](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))

4. Adapter and quality trimming ([`Trim Galore!`](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/))

5. Alignment to both target and spike-in genomes ([`Bowtie 2`](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml))

6. Filter on quality, sort and index alignments ([`SAMtools`](https://sourceforge.net/projects/samtools/files/samtools/))

7. Duplicate read marking ([`picard MarkDuplicates`](https://broadinstitute.github.io/picard/))

8. Create bedGraph files ([`BEDTools`](https://github.com/arq5x/bedtools2/)

9. Create bigWig coverage files ([`bedGraphToBigWig`](http://hgdownload.soe.ucsc.edu/admin/exe/))

10. Peak calling specifically tailored for low background noise ([`SEACR`](https://github.com/FredHutch/SEACR))

11. Quality control and analysis:

1. Merge re-sequenced FastQ files ([`cat`](http://www.linfo.org/cat.html))

2. Read QC ([`FastQC`](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))

3. Adapter and quality trimming ([`Trim Galore!`](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/))

4. Alignment to both target and spike-in genomes ([`Bowtie 2`](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml))

5. Filter on quality, sort and index alignments ([`SAMtools`](https://sourceforge.net/projects/samtools/files/samtools/))

6. Duplicate read marking ([`picard MarkDuplicates`](https://broadinstitute.github.io/picard/))

7. Create bedGraph files ([`BEDTools`](https://github.com/arq5x/bedtools2/)

8. Create bigWig coverage files ([`bedGraphToBigWig`](http://hgdownload.soe.ucsc.edu/admin/exe/))

9. Peak calling specifically tailored for low background noise ([`SEACR`](https://github.com/FredHutch/SEACR))

10. Quality control and analysis:

    1. Alignment, fragment length and peak analysis and replicate reproducibility ([`Python`](https://www.python.org/))

    2. Differential peak analysis ([`DESeq2`](https://bioconductor.org/packages/release/bioc/html/DESeq2.html))

    3. Heatmap peak analysis ([`deepTools`](https://github.com/deeptools/deepTools/))

12. Genome browser session ([`IGV`](https://software.broadinstitute.org/software/igv/))

13. Present QC for raw read, alignment and duplicate reads ([`MultiQC`](http://multiqc.info/))

11. Genome browser session ([`IGV`](https://software.broadinstitute.org/software/igv/))

12. Present QC for raw read, alignment and duplicate reads ([`MultiQC`](http://multiqc.info/))

## Quick Start

            --genome GRCh37

        ```

    * Typical command for downloading public data:

        ```bash

        nextflow run nf-core/cutandrun \

            --public_data_ids ids.txt \

            -profile <docker/singularity/podman/conda/institute>

        ```

    > **NB:** The commands to obtain public data and to run the main arm of the pipeline are completely independent. This is intentional because it allows you to download all of the raw data in an initial pipeline run (`results/public_data/`) and then to curate the auto-created samplesheet based on the available sample metadata before you run the pipeline again properly.

See [usage docs](https://nf-co.re/cutandrun/usage) for all of the available options when running the pipeline.

## Documentation

#!/usr/bin/env python

import os

import re

import sys

import csv

import errno

import requests

import argparse

## Example ids supported by this script

SRA_IDS = ['PRJNA63463', 'SAMN00765663', 'SRA023522', 'SRP003255', 'SRR390278', 'SRS282569', 'SRX111814']

ENA_IDS = ['ERA2421642', 'ERP120836', 'ERR674736', 'ERS4399631', 'ERX629702', 'PRJEB7743', 'SAMEA3121481']

GEO_IDS = ['GSE18729', 'GSM465244']

ID_REGEX = r'^[A-Z]+'

PREFIX_LIST = sorted(list(set([re.search(ID_REGEX,x).group() for x in SRA_IDS + ENA_IDS + GEO_IDS])))

def parse_args(args=None):

    Description = 'Download and create a run information metadata file from SRA/ENA/GEO identifiers.'

    Epilog = 'Example usage: python fetch_sra_runinfo.py <FILE_IN> <FILE_OUT>'

    parser = argparse.ArgumentParser(description=Description, epilog=Epilog)

    parser.add_argument('FILE_IN', help="File containing database identifiers, one per line.")

    parser.add_argument('FILE_OUT', help="Output file in tab-delimited format.")

    parser.add_argument('-pl', '--platform', type=str, dest="PLATFORM", default='', help="Comma-separated list of platforms to use for filtering. Accepted values = 'ILLUMINA', 'OXFORD_NANOPORE' (default: '').")

    parser.add_argument('-ll', '--library_layout', type=str, dest="LIBRARY_LAYOUT", default='', help="Comma-separated list of library layouts to use for filtering. Accepted values = 'SINGLE', 'PAIRED' (default: '').")

    return parser.parse_args(args)

def validate_csv_param(param,valid_vals,param_desc):

    valid_list = []

    if param:

        user_vals = param.split(',')

        intersect = list(set(user_vals) & set(valid_vals))

        if len(intersect) == len(user_vals):

            valid_list = intersect

        else:

            print("ERROR: Please provide a valid {} parameter!\nProvided values = {}\nAccepted values = {}".format(param_desc,param,','.join(validVals)))

            sys.exit(1)

    return valid_list

def make_dir(path):

    if not len(path) == 0:

        try:

            os.makedirs(path)

        except OSError as exception:

            if exception.errno != errno.EEXIST:

                raise

def fetch_url(url,encoding='utf-8'):

    try:

        r = requests.get(url)

    except requests.exceptions.RequestException as e:

        raise SystemExit(e)

    if r.status_code != 200:

        print("ERROR: Connection failed\nError code '{}'".format(r.status_code))

        sys.exit(1)

    return r.content.decode(encoding).splitlines()

def id_to_srx(db_id):

    ids = []

    url = 'https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term={}'.format(db_id)

    for row in csv.DictReader(fetch_url(url), delimiter=','):

        ids.append(row['Experiment'])

    return ids

def id_to_erx(db_id):

    ids = []

    fields = ['run_accession', 'experiment_accession']

    url = 'http://www.ebi.ac.uk/ena/data/warehouse/filereport?accession={}&result=read_run&fields={}'.format(db_id,','.join(fields))

    for row in csv.DictReader(fetch_url(url), delimiter='\t'):

        ids.append(row['experiment_accession'])

    return ids

def gse_to_srx(db_id):

    ids = []

    url = 'https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc={}&targ=gsm&view=data&form=text'.format(db_id)

    gsm_ids = [x.split('=')[1].strip() for x in fetch_url(url) if x.find('GSM') != -1]

    for gsm_id in gsm_ids:

        ids += id_to_srx(gsm_id)

    return ids

def get_ena_fields():

    fields = []

    url = 'https://www.ebi.ac.uk/ena/portal/api/returnFields?dataPortal=ena&format=tsv&result=read_run'

    for row in csv.DictReader(fetch_url(url), delimiter='\t'):

        fields.append(row['columnId'])

    return fields

def fetch_sra_runinfo(file_in,file_out,platform_list=[],library_layout_list=[]):

    total_out = 0

    seen_ids = []; run_ids = []

    header = []

    make_dir(os.path.dirname(file_out))

    ena_fields = get_ena_fields()

    with open(file_in,"r") as fin, open(file_out,"w") as fout:

        for line in fin:

            db_id = line.strip()

            match = re.search(ID_REGEX, db_id)

            if match:

                prefix = match.group()

                if prefix in PREFIX_LIST:

                    if not db_id in seen_ids:

                        ids = [db_id]

                        ## Resolve/expand these ids against GEO URL

                        if prefix in ['GSE']:

                            ids = gse_to_srx(db_id)

                        ## Resolve/expand these ids against SRA URL

                        elif prefix in ['GSM', 'PRJNA', 'SAMN', 'SRR']:

                            ids = id_to_srx(db_id)

                        ## Resolve/expand these ids against ENA URL

                        elif prefix in ['ERR']:

                            ids = id_to_erx(db_id)

                        ## Resolve/expand to get run identifier from ENA and write to file

                        for id in ids:

                            url = 'http://www.ebi.ac.uk/ena/data/warehouse/filereport?accession={}&result=read_run&fields={}'.format(id,','.join(ena_fields))

                            csv_dict = csv.DictReader(fetch_url(url), delimiter='\t')

                            for row in csv_dict:

                                run_id = row['run_accession']

                                if not run_id in run_ids:

                                    write_id = True

                                    if platform_list:

                                        if row['instrument_platform'] not in platform_list:

                                            write_id = False

                                    if library_layout_list:

                                        if row['library_layout'] not in library_layout_list:

                                            write_id = False

                                    if write_id:

                                        if total_out == 0:

                                            header = sorted(row.keys())

                                            fout.write('{}\n'.format('\t'.join(sorted(header))))

                                        else:

                                            if header != sorted(row.keys()):

                                                print("ERROR: Metadata columns do not match for id {}!\nLine: '{}'".format(run_id,line.strip()))

                                                sys.exit(1)

                                        fout.write('{}\n'.format('\t'.join([row[x] for x in header])))

                                        total_out += 1

                                    run_ids.append(run_id)

                        seen_ids.append(db_id)

                        if not ids:

                            print("ERROR: No matches found for database id {}!\nLine: '{}'".format(db_id,line.strip()))

                            sys.exit(1)

                else:

                    id_str = ', '.join([x + "*" for x in PREFIX_LIST])

                    print("ERROR: Please provide a valid database id starting with {}!\nLine: '{}'".format(id_str,line.strip()))

                    sys.exit(1)

            else:

                id_str = ', '.join([x + "*" for x in PREFIX_LIST])

                print("ERROR: Please provide a valid database id starting with {}!\nLine: '{}'".format(id_str,line.strip()))

                sys.exit(1)

def main(args=None):

    args = parse_args(args)

    platform_list = validate_csv_param(args.PLATFORM,valid_vals=['ILLUMINA'],param_desc='--platform')

    library_layout_list = validate_csv_param(args.LIBRARY_LAYOUT,valid_vals=['SINGLE', 'PAIRED'],param_desc='--library_layout')

    fetch_sra_runinfo(args.FILE_IN,args.FILE_OUT,platform_list,library_layout_list)

if __name__ == '__main__':

    sys.exit(main())

#!/usr/bin/env python

import os

import sys

import errno

import argparse

import collections

def parse_args(args=None):

    Description = "Create samplesheet with FTP download links and md5ums from sample information obtained via 'sra_ids_to_runinfo.py' script."

    Epilog = 'Example usage: python sra_runinfo_to_ftp.py <FILES_IN> <FILE_OUT>'

    parser = argparse.ArgumentParser(description=Description, epilog=Epilog)

    parser.add_argument('FILES_IN', help="Comma-separated list of metadata file created from 'sra_ids_to_runinfo.py' script.")

    parser.add_argument('FILE_OUT', help="Output file containing paths to download FastQ files along with their associated md5sums.")

    return parser.parse_args(args)

def make_dir(path):

    if not len(path) == 0:

        try:

            os.makedirs(path)

        except OSError as exception:

            if exception.errno != errno.EEXIST:

                raise

def parse_sra_runinfo(file_in):

    runinfo_dict = {}

    with open(file_in, "r") as fin:

        header = fin.readline().strip().split('\t')

        for line in fin:

            line_dict   = dict(zip(header,line.strip().split('\t')))

            line_dict   = collections.OrderedDict(sorted(list(line_dict.items())))

            run_id      = line_dict['run_accession']

            exp_id      = line_dict['experiment_accession']

            library     = line_dict['library_layout']

            fastq_files = line_dict['fastq_ftp']

            fastq_md5   = line_dict['fastq_md5']

            print(line_dict)

            db_id = exp_id

            sample_dict = collections.OrderedDict()

            if library == 'SINGLE':

                sample_dict = collections.OrderedDict([('fastq_1',''), ('fastq_2',''), ('md5_1',''), ('md5_2',''), ('single_end','true')])

                if fastq_files:

                    sample_dict['fastq_1']  = fastq_files

                    sample_dict['md5_1']    = fastq_md5

                else:

                    ## In some instances FTP links don't exist for FastQ files

                    ## These have to be downloaded via fastq-dump / fasterq-dump / parallel-fastq-dump via the run id

                    db_id = run_id

            elif library == 'PAIRED':

                sample_dict = collections.OrderedDict([('fastq_1',''), ('fastq_2',''), ('md5_1',''), ('md5_2',''), ('single_end','false')])

                if fastq_files:

                    fq_files = fastq_files.split(';')[-2:]

                    fq_md5   = fastq_md5.split(';')[-2:]

                    if len(fq_files) == 2:

                        if fq_files[0].find('_1.fastq.gz') != -1 and fq_files[1].find('_2.fastq.gz') != -1:

                            sample_dict['fastq_1'] = fq_files[0]

                            sample_dict['fastq_2'] = fq_files[1]

                            sample_dict['md5_1']   = fq_md5[0]

                            sample_dict['md5_2']   = fq_md5[1]

                        else:

                            print("Invalid FastQ files found for database id:'{}'!.".format(run_id))

                    else:

                        print("Invalid number of FastQ files ({}) found for paired-end database id:'{}'!.".format(len(fq_files), run_id))

                else:

                    db_id = run_id

            if sample_dict:

                sample_dict.update(line_dict)

                if db_id not in runinfo_dict:

                    runinfo_dict[db_id] = [sample_dict]

                else:

                    if sample_dict in runinfo_dict[db_id]:

                        print("Input run info file contains duplicate rows!\nLine: '{}'".format(line))

                    else:

                        runinfo_dict[db_id].append(sample_dict)

    return runinfo_dict

def sra_runinfo_to_ftp(files_in,file_out):

    samplesheet_dict = {}

    for file_in in files_in:

        runinfo_dict = parse_sra_runinfo(file_in)

        for db_id in runinfo_dict.keys():

            if db_id not in samplesheet_dict:

                samplesheet_dict[db_id] = runinfo_dict[db_id]

            else:

                print("Duplicate sample identifier found!\nID: '{}'".format(db_id))

    ## Write samplesheet with paths to FastQ files and md5 sums

    if samplesheet_dict:

        out_dir = os.path.dirname(file_out)

        make_dir(out_dir)

        with open(file_out, "w") as fout:

            header = ['id'] + list(samplesheet_dict[list(samplesheet_dict.keys())[0]][0].keys())

            fout.write("\t".join(header) + "\n")

            for db_id in sorted(samplesheet_dict.keys()):

                for idx,val in enumerate(samplesheet_dict[db_id]):

                    fout.write('\t'.join(["{}_T{}".format(db_id,idx+1)] + [val[x] for x in header[1:]]) + '\n')

def main(args=None):

    args = parse_args(args)

    sra_runinfo_to_ftp([x.strip() for x in args.FILES_IN.split(',')], args.FILE_OUT)

if __name__ == '__main__':

    sys.exit(main())

params {

    modules {

        "sra_ids_to_runinfo" {

            publish_dir     = "public_data"

            publish_files   = ["tsv":"runinfo"]

        }

        "sra_runinfo_to_ftp" {

            publish_dir     = "public_data"

            publish_files   = ["tsv":"runinfo"]

        }

        "sra_fastq_ftp" {

            publish_dir     = "public_data"

            publish_files   = ["fastq.gz":"", "md5":"md5"]

            args            = "-C - --max-time 1200"

        }

        "sra_to_samplesheet" {

            publish_dir     = "public_data"

            publish_files   = false

        }

        "sra_merge_samplesheet" {

            publish_dir     = "public_data"

        }

        "bowtie2_index" {

            publish_dir   = "genome/index"

        }

/*

 * -------------------------------------------------

 *  Nextflow config file for running tests

 * -------------------------------------------------

 * Defines bundled input files and everything required

 * to run a fast and simple test. Use as follows:

 *   nextflow run nf-core/cutandrun -profile test_sra,<docker/singularity>

 */

params {

  config_profile_name        = 'Public data download test profile'

  config_profile_description = 'Minimal test dataset to check pipeline function when downloading data via the ENA'

  // Limit resources so that this can run CI

  max_cpus   = 2

  max_memory = 6.GB

  max_time   = 6.h

  // Input data

  public_data_ids = 'https://raw.githubusercontent.com/nf-core/test-datasets/rnaseq/samplesheet/public_database_ids.txt'

}

 No newline at end of file