snakemake 轉(zhuǎn)錄組上游pipeline

trim-galore + hisat2 + featureCounts

##工作的路徑別名，主要作用引用別名浴讯，省代碼帮非，同時(shí)更加簡潔
base_path = "/path/"
##所有樣本名的集合司蔬，用于后面的批量調(diào)用
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"

# 構(gòu)建索引
# 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>: 設(shè)定Phred quality score閾值，默認(rèn)為20棚愤。
# --phred33: 選擇-phred33或者-phred64，表示測序平臺(tái)使用的Phred quality score。
# --adapter: 輸入adapter序列徘溢。也可以不輸入，Trim Galore!會(huì)自動(dòng)尋找可能性最高的平臺(tái)對應(yīng)的adapter捆探。自動(dòng)搜選的平臺(tái)三個(gè)然爆，
# 也直接顯式輸入這三種平臺(tái)，即--illumina黍图、--nextera和--small_rna曾雕。
# --stringency <int>: 設(shè)定可以忍受的前后adapter重疊的堿基數(shù)，默認(rèn)為1（非持唬苛刻）剖张。可以適度放寬揩环，因?yàn)楹笠粋€(gè)adapter幾乎不可能被測序儀讀到搔弄。
# --length <int>: 設(shè)定輸出reads長度閾值，小于設(shè)定值會(huì)被拋棄丰滑。
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>: 線程數(shù)目
# -x <string>: 參考基因組索引的basename顾犹，即前綴名
# -1 <string>: 雙端測序的read1 list ，若為list褒墨，使用逗號(hào)隔開炫刷，名字與2要匹配，如-1 flyA_1.fq,flyB_1.fq
# -2 <string>: 雙端測序的read2 list 郁妈，若為list柬唯，使用逗號(hào)隔開，名字與1要匹配圃庭，如-2 flyA_2.fq,flyB_2.fq
# -S <string>: SAM寫入的文件名锄奢，默認(rèn)寫入到標(biāo)準(zhǔn)輸出中
# --dta: 注意Ｊ纭！拘央！在下游使用stringtie組裝的時(shí)候一定要在hisat中設(shè)置這個(gè)參數(shù)Ｍ科ā！灰伟！
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: 該參數(shù)設(shè)置輸出 BAM 格式拆又，默認(rèn)下輸出是 SAM 格式文件
# -h print header for the output: 默認(rèn)下輸出的文件不帶 header，該參數(shù)設(shè)定輸出文件時(shí)帶 header 信息
# -S input is SAM: 默認(rèn)下輸入是 BAM 文件栏账，若是輸入是 SAM 文件帖族，則最好加該參數(shù)，否則有時(shí)候會(huì)報(bào)錯(cuò)挡爵。
# -@ Number of additional threads to use [0]: 指使用的線程數(shù)
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>: 指使用的線程數(shù)
# -o <string>: 輸出文件的名字竖般，輸出文件的內(nèi)容為read 的統(tǒng)計(jì)數(shù)目
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的情況中，會(huì)統(tǒng)計(jì)fragment而不統(tǒng)計(jì)read
# -T <int>: 線程數(shù)目茶鹃，1~32
# -t <string>: 設(shè)置feature-type涣雕，-t指定的必須是gtf中有的feature，同時(shí)read只有落到這些feature上才會(huì)被統(tǒng)計(jì)到闭翩，默認(rèn)是“exon”
# -g <string>: 當(dāng)參考的gtf提供的時(shí)候挣郭，我們需要提供一個(gè)id identifier 來將feature水平的統(tǒng)計(jì)匯總為meta-feature水平的統(tǒng)計(jì)，
# 默認(rèn)為gene_id疗韵，注意兑障！選擇gtf中提供的id identifier!!!

# -a <string>: 參考gtf文件名，支持Gzipped文件格式
# -o <string>: 輸出文件的名字蕉汪，輸出文件的內(nèi)容為read 的統(tǒng)計(jì)數(shù)目
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"