snakemake 转录组上游pipeline

trim-galore + hisat2 + featureCounts

##工作的路径别名，主要作用引用别名，省代码，同时更加简洁
base_path = "/path/"
##所有样本名的集合，用于后面的批量调用
tissue_type = {"1-ICM-P16_FRAS210019583-1r", "1-ICM-P17_FRAS210019584-2r", "1-ICM-P18_FRAS210019585-1r",
               "3-ICM-P16_FRAS210019580-1r", "3-ICM-P17_FRAS210019581-1r", "3-ICM-P18_FRAS210019582-2r",
               "IPS-4-P19_FRAS210019586-1r", "IPS-4-P20_FRAS210019587-2r", "IPS-4-P21_FRAS210019588-1r",
               "X-2-P5_FRAS210019577-1r", "X-2-P6_FRAS210019578-1r", "X-2-P7_FRAS210019579-1r", "OFC_FRAS210019589-2r"}

hisat2_index = "/path/to/hisat2_genome/genome"
annotation_gtf = "/path/to/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gff"

# 构建索引
# hisat2-build -p 20 genome.fa genome

rule all:
    input:
        expand("{base_path}/bam/{tissue_type}_hisat_sorted.bam",
            base_path=base_path, tissue_type=tissue_type),
        expand("{base_path}/bam/{tissue_type}_hisat_sorted.bam.bai",
            base_path=base_path, tissue_type=tissue_type),
        expand("{base_path}/count/counts.txt", base_path=base_path),


# --quality <int>: 设定Phred quality score阈值，默认为20。
# --phred33: 选择-phred33或者-phred64，表示测序平台使用的Phred quality score。
# --adapter: 输入adapter序列。也可以不输入，Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，
# 也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
# --stringency <int>: 设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
# --length <int>: 设定输出reads长度阈值，小于设定值会被抛弃。
rule trim_galore:
    input:
        "{base_path}/raw/{tissue_type}_{rep_num}_1.fq.gz",
        "{base_path}/raw/{tissue_type}_{rep_num}_2.fq.gz"
    output:
        "{base_path}/clean/{tissue_type}_{rep_num}_1.clean.fq.gz",
        "{base_path}/clean/{tissue_type}_{rep_num}_2.clean.fq.gz"
    log:
        "{base_path}/clean/{tissue_type}_{rep_num}.log"
    shell:
        "trim_galore -q 20 --phred33 --length 50 -e 0.1 --stringency 3 \
        -o {output[0]} {output[1]} {input[0]} {input[1]} > {log} 2>&1"


# -p <int>: 线程数目
# -x <string>: 参考基因组索引的basename，即前缀名
# -1 <string>: 双端测序的read1 list ，若为list，使用逗号隔开，名字与2要匹配，如-1 flyA_1.fq,flyB_1.fq
# -2 <string>: 双端测序的read2 list ，若为list，使用逗号隔开，名字与1要匹配，如-2 flyA_2.fq,flyB_2.fq
# -S <string>: SAM写入的文件名，默认写入到标准输出中
# --dta: 注意！！！在下游使用stringtie组装的时候一定要在hisat中设置这个参数！！！
rule hisat2:
    input:
        "{base_path}/clean/{tissue_type}_1.clean.fq.gz",
        "{base_path}/clean/{tissue_type}_2.clean.fq.gz"
    output:
        temp("{base_path}/bam/{tissue_type}_hisat.sam")
    log:
        "{base_path}/bam/{tissue_type}_hisat.log"
    threads:
        8
    shell:
        "hisat2 -p 8 -x {hisat2_index} -1 {input[0]} -2 {input[1]} -S {output[0]} 1>{log} 2>&1"


# -b output BAM: 该参数设置输出 BAM 格式，默认下输出是 SAM 格式文件
# -h print header for the output: 默认下输出的文件不带 header，该参数设定输出文件时带 header 信息
# -S input is SAM: 默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
# -@ Number of additional threads to use [0]: 指使用的线程数
rule sam2bam:
    input:
        "{base_path}/bam/{tissue_type}_hisat.sam"
    output:
        temp("{base_path}/bam/{tissue_type}_hisat.bam")
    threads:
        8
    shell:
        "samtools view -b -S -h -@ 8 {input[0]} > {output[0]}"


# -@ <int>: 指使用的线程数
# -o <string>: 输出文件的名字，输出文件的内容为read 的统计数目
rule bam_sort:
    input:
        "{base_path}/bam/{tissue_type}_hisat.bam"
    output:
        "{base_path}/bam/{tissue_type}_hisat_sorted.bam"
    log:
        "{base_path}/bam/{tissue_type}_bam_sort.log"
    threads:
        8
    shell:
        "samtools sort -@ 8 -o {output[0]} {input[0]} 1>{log} 2>&1"


rule bam_index:
    input:
        "{base_path}/bam/{tissue_type}_hisat_sorted.bam"
    output:
        "{base_path}/bam/{tissue_type}_hisat_sorted.bam.bai"
    shell:
        "samtools index {input[0]} {output[0]}"


# -p: 只能用在paired-end的情况中，会统计fragment而不统计read
# -T <int>: 线程数目，1~32
# -t <string>: 设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”
# -g <string>: 当参考的gtf提供的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，
# 默认为gene_id，注意！选择gtf中提供的id identifier!!!

# -a <string>: 参考gtf文件名，支持Gzipped文件格式
# -o <string>: 输出文件的名字，输出文件的内容为read 的统计数目
rule featureCounts:
    input:
        expand("{base_path}/bam/{tissue_type}_hisat_sorted.bam", base_path=base_path,
            tissue_type=tissue_type)
    output:
        "{base_path}/count/counts.txt"
    log:
        "{base_path}/count/featureCounts.log"
    threads:
        8
    shell:
        "ls {base_path}/bam/*_hisat_sorted.bam | xargs featureCounts -a {annotation_gtf} -o {output[0]} \
        -p -T 8 -F GFF -t gene -g ID 1>{log} 2>&1"