Snakefile 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704
  1. import os
  2. import yaml
  3. configfile: "config.yaml"
  4. ##########################################################
  5. # 全局变量和样本信息
  6. ##########################################################
  7. ## 软件主目录
  8. PEAKSNAKE_HOME = config["PEAKSNAKE_HOME"] if config["PEAKSNAKE_HOME"] else os.getcwd()
  9. PEAK_TYPE = config["peak_type"]
  10. SEQ_TYPE = config["seq_type"]
  11. PEAK_SELECTION = config["peak_selection"]
  12. GENOME = config["genome"]
  13. GTF = config["gtf"]
  14. ## 样本信息变量
  15. # 字典:REPLICATE to INPUT
  16. REPLICATE_TO_INPUT = {k: config['samples'][k] for k in sorted(config['samples'])}
  17. # 列表:所有 REPLICATES
  18. REPLICATES = sorted(list(set(REPLICATE_TO_INPUT.keys())))
  19. # 列表:所有 INPUTS
  20. INPUTS = sorted(list(set(REPLICATE_TO_INPUT.values())))
  21. INPUTS = [] if INPUTS == [None] else INPUTS
  22. # 字典: SAMPLE 与 rep1 rep2 rep2 对应关系
  23. SAMPLE_TO_REPLICATE = {}
  24. for s in REPLICATES:
  25. name, rep = '_'.join(s.split('_')[:-1]), s.split('_')[-1]
  26. SAMPLE_TO_REPLICATE.setdefault(name, []).append(rep)
  27. SAMPLE_TO_REPLICATE = {k: sorted(v) for k, v in SAMPLE_TO_REPLICATE.items()}
  28. ## 生成样本信息表
  29. with open("sample_sheet.csv", 'w') as f:
  30. f.write("SampleID,ControlID,Tissue,Factor,Condition,Treatment,Replicate,bamReads,Peaks,bamControl,PeakCaller\n")
  31. for sample, control in REPLICATE_TO_INPUT.items():
  32. sample_parts = sample.split('_')
  33. factor = sample_parts[0] # 提取样本ID中的Factor
  34. tissue = "NA"
  35. treatment = sample_parts[1] # 在这个例子中,tissue和condition是相同的
  36. condition = factor + "_" + treatment
  37. replicate = sample_parts[2].replace("rep", "") # 将"rep1"和"rep2"转换为"1"和"2"
  38. if control:
  39. control_parts = control.split('_')
  40. control_id = "_".join(control_parts[:2]) # 构建 ControlID
  41. bamControl = f"clean_bams/{control}_final.bam"
  42. else:
  43. control_id = "NA"
  44. bamControl = "NA"
  45. bamReads = f"clean_bams/{sample}_final.bam"
  46. peaks = f"clean_peaks/cutoff/{sample}_peaks.{PEAK_TYPE}Peak"
  47. f.write(f"{sample},{control_id},{tissue},{factor},{condition},{treatment},{replicate},{bamReads},{peaks},{bamControl},bed\n")
  48. ##########################################################
  49. # rule all: 最终想要生成的文件
  50. ##########################################################
  51. rule all:
  52. input:
  53. #####################################
  54. # 从 fastq 到 peaks
  55. #####################################
  56. # 最终比对结果
  57. expand("clean_bams/{replicate}_final.bam", replicate = REPLICATES + INPUTS),
  58. # bw 文件, deeptools 可视化
  59. expand("clean_bams/{replicate}.bw", replicate = REPLICATES + INPUTS),
  60. "deeptools/sample_correlation.pdf" if len(REPLICATES) > 2 else [],
  61. "deeptools/tss_heatmap.pdf",
  62. "deeptools/tss_tes_heatmap.pdf",
  63. # raw peak 结果
  64. expand("raw_peaks/{replicate}_peaks.{PEAK_TYPE}Peak", replicate = REPLICATES, PEAK_TYPE = PEAK_TYPE),
  65. # cutoff analysis
  66. expand("raw_peaks/{replicate}_cutoff_analysis.pdf", replicate = REPLICATES),
  67. # clean peak 结果
  68. expand("clean_peaks/cutoff/{replicate}_peaks.{PEAK_TYPE}Peak", replicate = REPLICATES, PEAK_TYPE = PEAK_TYPE),
  69. # 通过 intersect or idr 进行 peak 筛选
  70. expand("clean_peaks/{m}/{sample}_peaks.{PEAK_TYPE}Peak", sample=SAMPLE_TO_REPLICATE.keys(), m = PEAK_SELECTION, PEAK_TYPE = PEAK_TYPE),
  71. # 提取序列
  72. expand("clean_peaks/{m}/{sample}_peaks.fa", sample=SAMPLE_TO_REPLICATE.keys(), m = PEAK_SELECTION),
  73. # 所有样本共识peaks
  74. "clean_peaks/merge/merged_peaks.bed" if len(SAMPLE_TO_REPLICATE) > 2 else [],
  75. #####################################
  76. # Peak 注释
  77. #####################################
  78. expand("peak_annotation/{sample}_peaks_allhits.txt", sample = SAMPLE_TO_REPLICATE.keys()) if config["uropa"]["do"] else [],
  79. expand("peak_annotation/{sample}_peaks_finalhits.txt", sample = SAMPLE_TO_REPLICATE.keys()) if config["uropa"]["do"] else [],
  80. #####################################
  81. # 定量分析
  82. #####################################
  83. "counts/merged_peaks.counts.matrix",
  84. "counts/merged_peaks.TMM.CPM.matrix",
  85. # 差异分析
  86. expand("diff_peaks/{contrast}.{m}.DE_results", contrast=config["contrasts"], m=config["diff_peaks"]["method"]) if config["diff_peaks"]["do"] and config["contrasts"] else [],
  87. expand("diff_peaks_annotation/{contrast}.bed6", contrast=config["contrasts"]) if config["diff_peaks"]["do"] and config["contrasts"] and config["uropa"]["do"] else [],
  88. #####################################
  89. # 质控结果
  90. #####################################
  91. # 测序数据质控表格
  92. "quality_control/fastp_summary.tsv",
  93. # 比对质控表格
  94. "quality_control/estimate_read_filtering.tsv",
  95. # cross correlation 质控表格
  96. "quality_control/cross_correlation_summary.tsv" if config["cross_correlation"]["do"] else [],
  97. ##########################################################
  98. # 使用fastp进行原始数据质控和过滤,自动判断 SE?PE?
  99. # raw_data => clean_data
  100. ##########################################################
  101. rule fastp_quality_control:
  102. input:
  103. fastq=["raw_data/{replicate}_R1.fastq.gz", "raw_data/{replicate}_R2.fastq.gz"] if SEQ_TYPE == 'PE' else ["raw_data/{replicate}_R1.fastq.gz"]
  104. output:
  105. fastq=["clean_data/{replicate}_R1.fastq.gz", "clean_data/{replicate}_R2.fastq.gz"] if SEQ_TYPE == 'PE' else ["clean_data/{replicate}_R1.fastq.gz"],
  106. html_report="clean_data/{replicate}_fastp.html",
  107. json_report="clean_data/{replicate}_fastp.json"
  108. log: "logs/{replicate}_fastp_quality_control.log"
  109. threads: 3
  110. singularity:
  111. PEAKSNAKE_HOME + "/sifs/commonTools_20231218.sif"
  112. params:
  113. fastq="-i raw_data/{replicate}_R1.fastq.gz -o clean_data/{replicate}_R1.fastq.gz -I raw_data/{replicate}_R2.fastq.gz -O clean_data/{replicate}_R2.fastq.gz" if SEQ_TYPE == 'PE' else "-i raw_data/{replicate}_R1.fastq.gz -o clean_data/{replicate}_R1.fastq.gz",
  114. fastp=config["fastp"]
  115. shell:
  116. """
  117. fastp {params.fastq} --html {output.html_report} --json {output.json_report} --thread {threads} {params.fastp} 1>{log} 2>&1
  118. """
  119. ##########################################################
  120. # 根据基因组构建 bowtie2 index
  121. ##########################################################
  122. rule bowtie2_build:
  123. input:
  124. genome = GENOME
  125. output:
  126. index_prefix = "ref/genome.1.bt2"
  127. log: "logs/bowtie2_build.log"
  128. singularity:
  129. PEAKSNAKE_HOME + "/sifs/bowtie2_2.5.2.sif"
  130. threads: 6
  131. shell:
  132. "bowtie2-build --threads {threads} {input.genome} ref/genome 1>{log} 2>&1"
  133. ##########################################################
  134. # 使用 bowtie2 将测序数据比对到参考基因组
  135. # fastq => sam
  136. ##########################################################
  137. rule bowtie2_align:
  138. input:
  139. fastq=["clean_data/{replicate}_R1.fastq.gz", "clean_data/{replicate}_R2.fastq.gz"] if SEQ_TYPE == 'PE' else ["clean_data/{replicate}_R1.fastq.gz"],
  140. genome_index=config["bowtie2_index"] + ".1.bt2" if config["bowtie2_index"] else "ref/genome.1.bt2"
  141. output:
  142. sam="raw_bams/{replicate}.sam"
  143. log: "logs/{replicate}_bowtie2_align.log"
  144. singularity:
  145. PEAKSNAKE_HOME + "/sifs/bowtie2_2.5.2.sif"
  146. threads: 4
  147. params:
  148. genome_index=config["bowtie2_index"] if config["bowtie2_index"] else "ref/genome",
  149. fastq="-1 clean_data/{replicate}_R1.fastq.gz -2 clean_data/{replicate}_R2.fastq.gz" if SEQ_TYPE == 'PE' else "-U clean_data/{replicate}_R1.fastq.gz",
  150. bowtie2=config["bowtie2"] if config["bowtie2"] else ""
  151. shell:
  152. """
  153. bowtie2 -p {threads} -x {params.genome_index} {params.fastq} -S {output.sam} {params.bowtie2} 1>{log} 2>&1
  154. """
  155. rule sort_bam:
  156. input:
  157. sam="raw_bams/{replicate}.sam"
  158. output:
  159. sorted_bam="raw_bams/{replicate}_sorted.bam"
  160. log: "logs/{replicate}_sort_bam.log"
  161. singularity:
  162. PEAKSNAKE_HOME + "/sifs/commonTools_20231218.sif"
  163. threads: 4
  164. shell:
  165. """
  166. samtools sort -@ {threads} -o {output.sorted_bam} {input.sam} 1>{log} 2>&1
  167. samtools index {output.sorted_bam}
  168. """
  169. rule sieve_alignment:
  170. input:
  171. bam="raw_bams/{replicate}_sorted.bam"
  172. output:
  173. bam="clean_bams/{replicate}_final.bam"
  174. log: "logs/{replicate}_sieve_alignment.log"
  175. threads: 4
  176. singularity:
  177. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  178. params:
  179. peaksnake_home=PEAKSNAKE_HOME,
  180. minMappingQuality=config["alignmentSieve"]["minMappingQuality"],
  181. blacklist="--blackListFileName " + config["alignmentSieve"]["blacklist"] if config["alignmentSieve"]["blacklist"] else [],
  182. extra=config["alignmentSieve"]["extra"] if config["alignmentSieve"]["extra"] else ""
  183. shell:
  184. """
  185. alignmentSieve --numberOfProcessors {threads} --bam {input.bam} --outFile {output.bam} --filterMetrics {log} --ignoreDuplicates --minMappingQuality {params.minMappingQuality} --samFlagExclude 260 {params.blacklist} {params.extra} 1>{log} 2>&1
  186. samtools index {output.bam}
  187. """
  188. rule estimate_read_filtering:
  189. input:
  190. bams=expand("raw_bams/{replicate}_sorted.bam", replicate=REPLICATES + INPUTS)
  191. output:
  192. stat="quality_control/estimate_read_filtering.tsv"
  193. threads: 4
  194. singularity:
  195. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  196. log: "logs/estimate_read_filtering.log"
  197. params:
  198. peaksnake_home=PEAKSNAKE_HOME,
  199. minMappingQuality=config["alignmentSieve"]["minMappingQuality"],
  200. blacklist="--blackListFileName " + config["alignmentSieve"]["blacklist"] if config["alignmentSieve"]["blacklist"] else "",
  201. extra=config["alignmentSieve"]["extra"] if config["alignmentSieve"]["extra"] else "",
  202. sampleLabels=REPLICATES + INPUTS
  203. shell:
  204. """
  205. estimateReadFiltering \
  206. --numberOfProcessors {threads} \
  207. --bam {input.bams} \
  208. --sampleLabels {params.sampleLabels} \
  209. --outFile {output.stat} \
  210. --ignoreDuplicates \
  211. --minMappingQuality {params.minMappingQuality} \
  212. --samFlagExclude 260 \
  213. --blackListFileName ref/blacklist.bed \
  214. {params.extra} 1>{log} 2>&1
  215. """
  216. ##########################################################
  217. # DeepTools 绘图
  218. ##########################################################
  219. rule convert_bam_to_bigwig:
  220. input:
  221. bam="clean_bams/{replicate}_final.bam"
  222. output:
  223. bw="clean_bams/{replicate}.bw"
  224. log: "logs/{replicate}_convert_bam_to_bigwig.log"
  225. singularity:
  226. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  227. threads: 2
  228. params:
  229. gsize=config["gsize"],
  230. bamCoverage = config["bamCoverage"] if config["bamCoverage"] else ""
  231. shell:
  232. """
  233. bamCoverage -p {threads} --bam {input.bam} -o {output.bw} --effectiveGenomeSize {params.gsize} {params.bamCoverage} 1>{log} 2>&1
  234. """
  235. rule summarize_multiple_bigwigs:
  236. input:
  237. bws=expand("clean_bams/{replicate}.bw", replicate=REPLICATES + INPUTS),
  238. tss_tes_bed=config["tss_tes_bed"]
  239. output:
  240. "clean_bams/bw_summary.gz"
  241. params:
  242. labels=REPLICATES + INPUTS
  243. log: "logs/summarize_multiple_bigwigs.log"
  244. singularity:
  245. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  246. threads: 10
  247. shell:
  248. """
  249. multiBigwigSummary BED-file \
  250. --bwfiles {input.bws} \
  251. --labels {params.labels} \
  252. --BED {input.tss_tes_bed} \
  253. -o {output} -p {threads} 1>{log} 2>&1
  254. """
  255. rule generate_correlation_plots:
  256. input:
  257. "clean_bams/bw_summary.gz"
  258. output:
  259. pdf="deeptools/sample_correlation.pdf",
  260. tab="deeptools/sample_correlation.tab"
  261. log: "logs/generate_correlation_plots.log"
  262. singularity:
  263. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  264. shell:
  265. """
  266. plotCorrelation -in {input} -o {output.pdf} \
  267. --corMethod pearson --whatToPlot heatmap \
  268. -min 0.5 \
  269. --plotTitle "Pearson Correlation of Samples" \
  270. --outFileCorMatrix {output.tab} 1>{log} 2>&1
  271. """
  272. rule plot_heatmap_reference_point:
  273. input:
  274. bws=expand("clean_bams/{replicate}.bw", replicate=REPLICATES + INPUTS),
  275. tss_tes_bed=config["tss_tes_bed"],
  276. tss_tes_shuffle_bed=config["tss_tes_shuffle_bed"]
  277. output:
  278. tss_matrix="deeptools/tss_matrix.gz",
  279. tss_heatmap="deeptools/tss_heatmap.pdf"
  280. params:
  281. labels=REPLICATES + INPUTS
  282. log: "logs/plot_heatmap_reference_point.log"
  283. singularity:
  284. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  285. threads: 10
  286. shell:
  287. """
  288. computeMatrix reference-point \
  289. -S {input.bws} \
  290. --samplesLabel {params.labels} \
  291. -R {input.tss_tes_bed} {input.tss_tes_shuffle_bed} \
  292. --referencePoint TSS \
  293. -b 5000 -a 5000 \
  294. --binSize 50 \
  295. -o {output.tss_matrix} \
  296. -p {threads} 1>{log} 2>&1
  297. plotHeatmap -m {output.tss_matrix} -o {output.tss_heatmap} --missingDataColor 0.5 1>>{log} 2>&1
  298. """
  299. rule plot_heatmap_scale_regions:
  300. input:
  301. bws=expand("clean_bams/{replicate}.bw", replicate=REPLICATES + INPUTS),
  302. tss_tes_bed=config["tss_tes_bed"],
  303. tss_tes_shuffle_bed=config["tss_tes_shuffle_bed"]
  304. output:
  305. tss_tes_matrix="deeptools/tss_tes_matrix.gz",
  306. tss_tes_heatmap="deeptools/tss_tes_heatmap.pdf"
  307. params:
  308. labels=REPLICATES + INPUTS
  309. log: "logs/plot_heatmap_scale_regions.log"
  310. singularity:
  311. PEAKSNAKE_HOME + "/sifs/deeptools_20231220.sif"
  312. threads: 10
  313. shell:
  314. """
  315. computeMatrix scale-regions \
  316. -S {input.bws} \
  317. --samplesLabel {params.labels} \
  318. -R {input.tss_tes_bed} {input.tss_tes_shuffle_bed} \
  319. --regionBodyLength 4000 \
  320. -b 2000 -a 2000 \
  321. --binSize 50 \
  322. -o {output.tss_tes_matrix} \
  323. -p {threads} 1>{log} 2>&1
  324. plotHeatmap -m {output.tss_tes_matrix} -o {output.tss_tes_heatmap} --missingDataColor 0.5 1>>{log} 2>&1
  325. """
  326. ##########################################################
  327. # 对每个 replicate:input call peak
  328. ##########################################################
  329. # 规则:MACS3 call peak
  330. rule callpeak_with_macs3:
  331. input:
  332. sorted_ip_bam="clean_bams/{replicate}_final.bam",
  333. sorted_input_bam=lambda wildcards: f"clean_bams/{config['samples'][wildcards.replicate]}_final.bam" if config['samples'][wildcards.replicate] else []
  334. output:
  335. Peak="raw_peaks/{replicate}_peaks." + PEAK_TYPE + "Peak",
  336. cutoff_analysis_txt="raw_peaks/{replicate}_cutoff_analysis.txt"
  337. log: "logs/{replicate}_callpeak_with_macs3.log"
  338. singularity:
  339. PEAKSNAKE_HOME + "/sifs/macs3_idr_20231218.sif"
  340. threads: 1
  341. params:
  342. control=lambda wildcards: f"-c clean_bams/{config['samples'][wildcards.replicate]}_final.bam" if config['samples'][wildcards.replicate] else "",
  343. format="BAMPE" if SEQ_TYPE == 'PE' else "BAM",
  344. gsize=config["gsize"],
  345. PEAK_TYPE="--broad" if PEAK_TYPE == "broad" else "",
  346. extra = config["macs3"]["extra"] if config["macs3"]["extra"] else ""
  347. shell:
  348. """
  349. macs3 callpeak -g {params.gsize} -t {input.sorted_ip_bam} {params.control} --name raw_peaks/{wildcards.replicate} --format {params.format} --keep-dup all --qvalue 0.075 --cutoff-analysis {params.PEAK_TYPE} {params.extra} 1>{log} 2>&1
  350. """
  351. rule plot_macs_cutoff_analysis:
  352. input:
  353. cutoff_analysis_txt="raw_peaks/{replicate}_cutoff_analysis.txt"
  354. output:
  355. cutoff_analysis_pdf="raw_peaks/{replicate}_cutoff_analysis.pdf"
  356. log: "logs/{replicate}_plot_macs_cutoff_analysis.log"
  357. shell:
  358. """
  359. python {PEAKSNAKE_HOME}/scripts/plot_macs_cutoff_analysis.py {input.cutoff_analysis_txt} {output} 1>{log} 2>&1
  360. """
  361. rule filter_peaks_by_qscore:
  362. input:
  363. Peak="raw_peaks/{replicate}_peaks." + PEAK_TYPE + "Peak"
  364. output:
  365. Peak="clean_peaks/cutoff/{replicate}_peaks." + PEAK_TYPE + "Peak"
  366. params:
  367. qscore=config["macs3"]["qscore"]
  368. shell:
  369. """
  370. awk '$9>{params.qscore}' {input.Peak} > {output.Peak}
  371. """
  372. rule select_peaks_norep:
  373. input:
  374. Peak="clean_peaks/cutoff/{sample}_rep1_peaks." + PEAK_TYPE + "Peak"
  375. output:
  376. Peak="clean_peaks/norep/{sample}_peaks." + PEAK_TYPE + "Peak"
  377. shell:
  378. """
  379. cp {input.Peak} {output.Peak}
  380. """
  381. rule select_peaks_by_intersect:
  382. input:
  383. Peak=lambda wildcards: expand("clean_peaks/cutoff/" + wildcards.sample + "_{rep}_peaks." + PEAK_TYPE + "Peak", rep=SAMPLE_TO_REPLICATE[wildcards.sample])
  384. output:
  385. Peak="clean_peaks/intersect/{sample}_peaks." + PEAK_TYPE + "Peak"
  386. params:
  387. min_overlap=config['intersect']['min_overlap']
  388. shell:
  389. """
  390. # 检查输入文件的数量
  391. num_peaks=$(echo {input.Peak} | wc -w)
  392. # 如果只有一个输入文件
  393. if [ "$num_peaks" -eq 1 ]; then
  394. cp {input.Peak[0]} {output.Peak}
  395. else
  396. # 复制第一个输入文件到临时文件
  397. cp {input.Peak[0]} {wildcards.sample}.temp.bed
  398. # 使用除第一个之外的所有输入文件
  399. echo {input.Peak} | tr ' ' '\\n' | awk 'NR>1' | while read bed; do
  400. bedtools intersect -f {params.min_overlap} -r -a {wildcards.sample}.temp.bed -b $bed -wa > {wildcards.sample}.temp2.bed
  401. mv {wildcards.sample}.temp2.bed {wildcards.sample}.temp.bed
  402. done
  403. # 创建一个中间的 peak list 文件
  404. cut -f 4 {wildcards.sample}.temp.bed > {wildcards.sample}.peak_lst
  405. # 使用中间的 peak list 文件和第一个输入文件生成最终输出
  406. awk 'NR==FNR {{a[$1]; next}} $4 in a' {wildcards.sample}.peak_lst {input.Peak[0]} > {output.Peak}
  407. # 清理临时文件
  408. rm {wildcards.sample}.temp.bed {wildcards.sample}.peak_lst
  409. fi
  410. """
  411. rule select_peaks_by_idr:
  412. input:
  413. rep1_peaks="raw_peaks/{sample}_rep1_peaks." + PEAK_TYPE + "Peak",
  414. rep2_peaks="raw_peaks/{sample}_rep2_peaks." + PEAK_TYPE + "Peak"
  415. output:
  416. true_rep_idr="clean_peaks/idr/{sample}_true_rep_idr.txt",
  417. idr_peaks="clean_peaks/idr/{sample}_peaks." + PEAK_TYPE + "Peak"
  418. log: "logs/{sample}_select_peaks_by_idr.log"
  419. singularity:
  420. PEAKSNAKE_HOME + "/sifs/macs3_idr_20231218.sif"
  421. threads: 1
  422. params:
  423. PEAK_TYPE=PEAK_TYPE,
  424. idr_threshold=config["idr"]["idr_threshold"]
  425. shell:
  426. """
  427. idr --samples {input.rep1_peaks} {input.rep2_peaks} --idr-threshold {params.idr_threshold} --output-file {output.true_rep_idr} --plot --input-file-type {PEAK_TYPE}Peak --rank p.value 1>{log} 2>&1
  428. awk -v OFS="\\t" 'BEGIN {{FS=OFS}} {{ $4=$1":"$2"_"$3; print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10 }}' {output.true_rep_idr} | sort -k1,1 -k2,2n -k3,3n > {output.idr_peaks}
  429. """
  430. ##########################################################
  431. # 提取序列
  432. ##########################################################
  433. rule extract_peak_sequence:
  434. input:
  435. Peak="clean_peaks/" + PEAK_SELECTION + "/{sample}_peaks." + PEAK_TYPE + "Peak"
  436. output:
  437. peak_fa="clean_peaks/" + PEAK_SELECTION + "/{sample}_peaks.fa"
  438. params:
  439. top_n=config["motif"]["top_n"],
  440. summit_flank=config["motif"]["summit_flank"],
  441. genome=config["genome"],
  442. summit_fa="clean_peaks/" + PEAK_SELECTION + "/{sample}_summit.fa"
  443. shell:
  444. """
  445. set +o pipefail;
  446. # Sorting and extracting fasta for peak
  447. sort -k8,8nr {input.Peak} | head -n {params.top_n} | bedtools getfasta -fi {params.genome} -bed - > {output.peak_fa}
  448. # Handling narrow peaks
  449. if [[ "{PEAK_TYPE}" == "narrow" ]]; then
  450. sort -k8,8nr {input.Peak} | head -n {params.top_n} | awk -v OFS='\\t' -v flank={params.summit_flank} '{{print $1, $2+$10-flank, $2+$10+flank+1}}' | awk '$2>=0' | bedtools getfasta -fi {params.genome} -bed - > {params.summit_fa}
  451. fi
  452. """
  453. ##########################################################
  454. # Peak 注释
  455. ##########################################################
  456. rule peak_annotation_with_uropa:
  457. input:
  458. Peak="clean_peaks/" + PEAK_SELECTION + "/{sample}_peaks." + PEAK_TYPE + "Peak",
  459. gtf=config["gtf"]
  460. output:
  461. allhits="peak_annotation/{sample}_peaks_allhits.txt",
  462. finalhits="peak_annotation/{sample}_peaks_finalhits.txt"
  463. log: "logs/{sample}_peak_annotation_with_uropa.log"
  464. params:
  465. feature=config["uropa"]["feature"],
  466. feature_anchor=config["uropa"]["feature_anchor"],
  467. distance=config["uropa"]["distance"],
  468. relative_location=config["uropa"]["relative_location"]
  469. singularity:
  470. PEAKSNAKE_HOME + "/sifs/uropa_4.0.3--pyhdfd78af_0.sif"
  471. shell:
  472. """
  473. uropa --bed {input.Peak} --gtf {input.gtf} --feature {params.feature} --feature-anchor {params.feature_anchor} --distance {params.distance} --relative-location {params.relative_location} --show-attributes gene_id --outdir peak_annotation 1>{log} 2>&1
  474. """
  475. ##########################################################
  476. # 生成所有 samples 共识 peaks
  477. ##########################################################
  478. rule merge_peaks:
  479. input:
  480. sample_peaks=expand("clean_peaks/" + PEAK_SELECTION + "/{sample}_peaks." + PEAK_TYPE + "Peak", sample=SAMPLE_TO_REPLICATE.keys())
  481. output:
  482. merged_peaks_bed="clean_peaks/merge/merged_peaks.bed",
  483. merged_peaks_saf="clean_peaks/merge/merged_peaks.saf"
  484. singularity:
  485. PEAKSNAKE_HOME + "/sifs/commonTools_20231218.sif"
  486. params:
  487. fai=config['genome'] + ".fai"
  488. shell:
  489. """
  490. mkdir -p clean_peaks/merge
  491. cat {input.sample_peaks} > clean_peaks/merge/cat_peaks.bed
  492. bedtools sort -i clean_peaks/merge/cat_peaks.bed -g {params.fai} > clean_peaks/merge/cat_sorted_peaks.bed
  493. bedtools merge -i clean_peaks/merge/cat_sorted_peaks.bed > {output.merged_peaks_bed}
  494. awk 'OFS="\t" {{print $1":"$2+1"-"$3, $1, $2+1, $3, "."}}' {output.merged_peaks_bed} > {output.merged_peaks_saf}
  495. """
  496. ##########################################################
  497. # 使用 feature counts 计算 replicate peak counts
  498. ##########################################################
  499. rule run_feature_counts:
  500. input:
  501. bam="clean_bams/{replicate}_final.bam",
  502. merged_peaks_saf="clean_peaks/merge/merged_peaks.saf"
  503. output:
  504. counts="counts/{replicate}.count",
  505. stat="counts/{replicate}.log"
  506. log: "logs/{replicate}_run_feature_counts.log"
  507. params:
  508. peaksnake_home=PEAKSNAKE_HOME,
  509. isPairedEnd="TRUE" if SEQ_TYPE == 'PE' else "FALSE"
  510. shell:
  511. """
  512. Rscript {params.peaksnake_home}/software/RunFeatureCounts/run-featurecounts.R -b {input.bam} --saf {input.merged_peaks_saf} --isPairedEnd {params.isPairedEnd} -o counts/{wildcards.replicate} 1>{log} 2>&1
  513. """
  514. ##########################################################
  515. # 合并生成 replicate counts 矩阵
  516. ##########################################################
  517. rule create_count_matrix:
  518. input:
  519. expand("counts/{replicate}.count", replicate=REPLICATES)
  520. output:
  521. counts_matrix="counts/merged_peaks.counts.matrix",
  522. cpm_matrix="counts/merged_peaks.TMM.CPM.matrix"
  523. log: "logs/create_count_matrix.log"
  524. params:
  525. peaksnake_home=PEAKSNAKE_HOME
  526. shell:
  527. """
  528. ls {input} >counts/count.fofn
  529. perl {params.peaksnake_home}/software/RunFeatureCounts/abundance_estimates_to_matrix.pl --est_method featureCounts --quant_files counts/count.fofn --out_prefix counts/merged_peaks 1>{log} 2>&1
  530. """
  531. ##########################################################
  532. # 使用 DESeq2/ edgeR 进行差异分析
  533. ##########################################################
  534. if config["contrasts"]:
  535. # 生成 samples.txt
  536. with open('samples.txt', 'w') as file:
  537. for ip in REPLICATES:
  538. # 提取 sample name
  539. sample_name = '_'.join(ip.split('_')[:2])
  540. file.write(f"{sample_name}\t{ip}\n")
  541. # 生成 contrasts.txt
  542. with open("contrasts.txt", "w") as file:
  543. file.write("\n".join([" ".join(item.split("_vs_")) for item in config["contrasts"]]))
  544. rule run_deseq2:
  545. input:
  546. "counts/merged_peaks.counts.matrix"
  547. output:
  548. expand("diff_peaks/{contrast}.DESeq2.DE_results", contrast=config["contrasts"])
  549. log: "logs/run_deseq2.log"
  550. singularity:
  551. PEAKSNAKE_HOME + "/sifs/trinity_20231218.sif"
  552. shell:
  553. """
  554. run_DE_analysis.pl -m {input} --method DESeq2 -s samples.txt --contrasts contrasts.txt -o diff_peaks 1>{log} 2>&1
  555. cd diff_peaks/
  556. rm *count_matrix *.MA_n_Volcano.pdf *.Rscript
  557. for file in merged_peaks.counts.matrix.*.DE_results; do
  558. echo -e "PeakID\tsampleA\tsampleB\tFoldChange\tPValue\tPadj" >"${{file#merged_peaks.counts.matrix.}}"
  559. grep -v "^sampleA" $file | awk 'BEGIN {{OFS="\t"}} {{ $7 = exp(log(2) * $7); print $1,$2,$3,$7,$10,$11}}' >> "${{file#merged_peaks.counts.matrix.}}"
  560. rm $file
  561. done
  562. cd ..
  563. """
  564. rule run_edgeR:
  565. input:
  566. "counts/merged_peaks.counts.matrix"
  567. output:
  568. expand("diff_peaks/{contrast}.edgeR.DE_results", contrast=config["contrasts"])
  569. log: "logs/run_edgeR.log"
  570. singularity:
  571. PEAKSNAKE_HOME + "/sifs/trinity_20231218.sif"
  572. shell:
  573. """
  574. run_DE_analysis.pl -m {input} --method edgeR -s samples.txt --contrasts contrasts.txt -o diff_peaks 1>{log} 2>&1
  575. cd diff_peaks/
  576. rm *count_matrix *.MA_n_Volcano.pdf *.Rscript
  577. for file in merged_peaks.counts.matrix.*.DE_results; do
  578. echo -e "PeakID\tsampleA\tsampleB\tFoldChange\tPValue\tPadj" >"${{file#merged_peaks.counts.matrix.}}"
  579. grep -v "^sampleA" $file | awk 'BEGIN {{OFS="\t"}} {{ $4 = exp(log(2) * $4); print $1,$2,$3,$4,$6,$7}}' >> "${{file#merged_peaks.counts.matrix.}}"
  580. rm $file
  581. done
  582. cd ..
  583. """
  584. rule diff_peaks_annotation_with_uropa:
  585. input:
  586. diff_peaks="diff_peaks/{contrast}.edgeR.DE_results",
  587. gtf=config["gtf"]
  588. output:
  589. diff_peaks_bed="diff_peaks_annotation/{contrast}.bed6",
  590. allhits="diff_peaks_annotation/{contrast}_allhits.txt",
  591. finalhits="diff_peaks_annotation/{contrast}_finalhits.txt"
  592. log: "logs/{contrast}_diff_peaks_annotation_with_uropa.log"
  593. params:
  594. diff_peaks_padj=config["diff_peaks"]["padj"],
  595. feature_anchor=config["uropa"]["feature_anchor"],
  596. distance=config["uropa"]["distance"],
  597. relative_location=config["uropa"]["relative_location"]
  598. singularity:
  599. PEAKSNAKE_HOME + "/sifs/uropa_4.0.3--pyhdfd78af_0.sif"
  600. shell:
  601. """
  602. awk 'BEGIN {{FS="\\t"; OFS="\t"}} NR>1 && $6 < {params.diff_peaks_padj} {{split($1, arr, "[:-]"); print arr[1], arr[2]-1, arr[3], $1, $4, "."}}' {input.diff_peaks} > {output.diff_peaks_bed}
  603. uropa --bed {output.diff_peaks_bed} --gtf {input.gtf} --feature gene --feature-anchor {params.feature_anchor} --distance {params.distance} --relative-location {params.relative_location} --outdir diff_peaks_annotation 1>{log} 2>&1
  604. """
  605. ##########################################################
  606. # 合并fastp质控结果生成表格
  607. ##########################################################
  608. rule combine_fastp_reports:
  609. input:
  610. expand("clean_data/{replicate}_fastp.json", replicate=REPLICATES + INPUTS)
  611. output:
  612. "quality_control/fastp_summary.tsv"
  613. log: "logs/combine_fastp_reports.log"
  614. params:
  615. peaksnake_home=PEAKSNAKE_HOME
  616. shell:
  617. """
  618. Rscript {params.peaksnake_home}/scripts/combine_fastp_report.R --input {input} --output quality_control/ 1>{log} 2>&1
  619. """
  620. ##########################################################
  621. # cross correlation
  622. ##########################################################
  623. rule cross_correlation:
  624. input:
  625. bam="clean_bams/{replicate}_final.bam"
  626. output:
  627. pdf="quality_control/cross_correlation/{replicate}.pdf",
  628. tsv="quality_control/cross_correlation/{replicate}.tsv"
  629. log: "logs/{replicate}_cross_correlation.log"
  630. singularity:
  631. PEAKSNAKE_HOME + "/sifs/phantompeakqualtools_20231218.sif"
  632. params:
  633. spp=config['cross_correlation']['spp'] if config['cross_correlation']['spp'] else ""
  634. threads: 4
  635. shell:
  636. """
  637. run_spp.R -p={threads} -c={input.bam} -savp={output.pdf} -out={output.tsv} -rf {params.spp} 1>{log} 2>&1
  638. """
  639. rule cross_correlation_summary:
  640. input:
  641. expand("quality_control/cross_correlation/{replicate}.tsv", replicate=REPLICATES + INPUTS)
  642. output:
  643. "quality_control/cross_correlation_summary.tsv"
  644. shell:
  645. """
  646. echo -e "Sample\\tTotalReadsFiltered\\tFragmentLength\\tCrossCorrelation\\tPhantomPeakLocation\\tPhantomPeakCorrelation\\tMinimumCrossCorrelationShift\\tMinimumCrossCorrelationValue\\tNSC\\tRSC\\tPhantomPeakQualityTag" >{output}
  647. cat {input} | sort -k 1 >>{output}
  648. """