# 準(zhǔn)備：
測(cè)序數(shù)據(jù)：fastq格式幌羞；RBFOX2 HepG2測(cè)序數(shù)據(jù) - RBFOX2
STAR的基因組索引：UCSC下載fasta格式
STAR的重復(fù)元件索引：RepBase下載fasta格式
barcodes文件：fasta格式
- PE: yeolabbarcodes_20170101.fasta
- SE: a_adapters.fasta 或者
InvRNA*_adapters.fasta
chrom.sizes文件：UCSC下載，tab分割硼一，第一列染色體名扯夭，第二列染色體長(zhǎng)度鳍贾，hg19 chrom.sizes
# 流程概述：
eclipdemux用于PE測(cè)序的樣本拆分和提取UMI；umi_tools用于SE測(cè)序的UMI提取
cutadapt修剪adapters
STAR比對(duì)到重復(fù)元件和篩選
比對(duì)篩選之后的測(cè)序數(shù)據(jù)到基因組
去除PCR產(chǎn)生的重復(fù)：SE測(cè)序使用umi_tools交洗，常規(guī)工具可以使用barcodecollapsepe.py
(paired-end only) Merges multiple inline barcodes and filters R1 (uses only R2 for peak calling)
Calls enriched peak regions (peak clusters) with CLIPPER
Uses size-matched input sample to normalize and calculate fold-change enrichment within enriched peak regions with custom perl scripts (overlap_peakfi_with_bam_PE.pl, peakscompress.pl)
## 1. eclipdemux用于PE測(cè)序的樣本拆分和提取UMI；umi_tools用于SE測(cè)序的UMI提取
# SE 鑒定UMI
umi_tools extract \
--random-seed 1 \
--bc-pattern NNNNNNNNNN \
--log EXAMPLE_SE.rep1_clip.---.--.metrics \ --stdin file_R1.fastq.gz \
--stdout EXAMPLE_SE.rep1.umi.r1.fq
# PE 樣本拆分和鑒定UMI
eclipdemux \
--metrics EXAMPLE_PE.rep1_clip.---.--.metrics \ --expectedbarcodeida C01 \ --expectedbarcodeidb D8f \
--fastq_1 file_R1.fastq.gz \
--fastq_2 file_R2.fastq.gz \
--newname rep2_clip \
--dataset EXAMPLE_PE \
--barcodesfile yeolabbarcodes_20170101.fasta \ --length 5
## ## 2. Cutadapt 去除adapters
因?yàn)榭赡馨l(fā)生兩次adapter連接事件橡淑，所以進(jìn)行兩次Cutadapt剪切
Cutadapt round 1:
cutadapt \
-f fastq \
--match-read-wildcards \
--times 1 \
-e 0.1 \
-O 1 \
--quality-cutoff 6 \
-m 18 \
-a NNNNNAGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \ -g CTTCCGATCTACAAGTT \
-g CTTCCGATCTTGGTCCT \
-A AACTTGTAGATCGGA \
-A AGGACCAAGATCGGA \
-A ACTTGTAGATCGGAA \
-A GGACCAAGATCGGAA \
-A CTTGTAGATCGGAAG \
-A GACCAAGATCGGAAG \
-A TTGTAGATCGGAAGA \
-A ACCAAGATCGGAAGA \
-A TGTAGATCGGAAGAG \
-A CCAAGATCGGAAGAG \
-A GTAGATCGGAAGAGC \
-A CAAGATCGGAAGAGC \
-A TAGATCGGAAGAGCG \
-A AAGATCGGAAGAGCG \
-A AGATCGGAAGAGCGT \
-A GATCGGAAGAGCGTC \
-A ATCGGAAGAGCGTCG \
-A TCGGAAGAGCGTCGT \
-A CGGAAGAGCGTCGTG \
-A GGAAGAGCGTCGTGT \
-o EXAMPLE_PE.rep2_clip.C01.r1.fqTr.fq \
-p EXAMPLE_PE.rep2_clip.C01.r2.fqTr.fq \ EXAMPLE_PE.rep2_clip.C01.r1.fq.gz \ EXAMPLE_PE.rep2_clip.C01.r2.fq.gz
Fastqc round 1
fastqc -t 2 --extract -k 7 EXAMPLE_PE.rep2_clip.C01.r1.fqTr.fq -o . fastqc -t 2 --extract -k 7 EXAMPLE_PE.rep2_clip.C01.r2.fqTr.fq –o .
Cutadapt round 2
cutadapt \
-f fastq \
--match-read-wildcards \
--times 1 \
-e 0.1 \
-O 5 \
--quality-cutoff 6 \
-m 18 \
-A AACTTGTAGATCGGA \
-A AGGACCAAGATCGGA \
-A ACTTGTAGATCGGAA \
-A GGACCAAGATCGGAA \
-A CTTGTAGATCGGAAG \
-A GACCAAGATCGGAAG \
-A TTGTAGATCGGAAGA \
-A ACCAAGATCGGAAGA \
-A TGTAGATCGGAAGAG \
-A CCAAGATCGGAAGAG \
-A GTAGATCGGAAGAGC \
-A CAAGATCGGAAGAGC \
-A TAGATCGGAAGAGCG \
-A AAGATCGGAAGAGCG \
-A AGATCGGAAGAGCGT \
-A GATCGGAAGAGCGTC \
-A ATCGGAAGAGCGTCG \
-A TCGGAAGAGCGTCGT \
-A CGGAAGAGCGTCGTG \
-A GGAAGAGCGTCGTGT \
-o EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.fq \ -p EXAMPLE_PE.rep2_clip.C01.r2.fqTrTr.fq \ EXAMPLE_PE.rep2_clip.C01.r1.fqTr.fq \ EXAMPLE_PE.rep2_clip.C01.r2.fqTr.fq
Fastqc round 2
fastqc -t 2 --extract -k 7 EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.fq -o . fastqc -t 2 --extract -k 7 EXAMPLE_PE.rep2_clip.C01.r2.fqTrTr.fq –o .
Fastq-sort
fastq-sort --id EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.fq > EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.sorted.fq
fastq-sort --id EXAMPLE_PE.rep2_clip.C01.r2.fqTrTr.fq > EXAMPLE_PE.rep2_clip.C01.r2.fqTrTr.sorted.fq
## 3. STAR比對(duì)到重復(fù)元件和篩選
#STAR rmRe
STAR \
--runMode alignReads \
--runThreadN 8 \
--genomeDir homo_sapiens_repbase_v2 \
--genomeLoad NoSharedMemory \
--alignEndsType EndToEnd \
--outSAMunmapped Within \
--outFilterMultimapNmax 30 \
--outFilterMultimapScoreRange 1 \
--outFileNamePrefix EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.sorted.STAR \ --outSAMtype BAM Unsorted \
--outFilterType BySJout \
--outBAMcompression 10 \
--outReadsUnmapped Fastx \
--outFilterScoreMin 10 \
--outSAMattrRGline ID:foo \
--outSAMattributes All \
--outSAMmode Full \
--outStd Log \
--readFilesIn EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.sorted.fq EXAMPLE_PE.rep2_clip.C01.r2.fqTrTr.sorted.fq
#Re-name files: re-name repeat-mapped outputs
mv EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.sorted.STARAligned.out.bam
EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-mapped.bam
mv EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.sorted.STARUnmapped.out.mate1 EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-unmapped.fq
mv EXAMPLE_PE.rep2_clip.C01.r1.fqTrTr.sorted.STARUnmapped.out.mate2 EXAMPLE_PE.rep2_clip.C01.r2.fq.repeat-unmapped.fq
# Fastq-sort
fastq-sort --id EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-unmapped.fq > EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-unmapped.sorted.fq fastq-sort --id EXAMPLE_PE.rep2_clip.C01.r2.fq.repeat-unmapped.fq > EXAMPLE_PE.rep2_clip.C01.r2.fq.repeat-unmapped.sorted.fq
## 4. 比對(duì)篩選之后的測(cè)序數(shù)據(jù)到基因組
#STAR genome mapping: Takes output from STAR rmRep. Maps unique reads to the human genome
STAR \
--runMode alignReads \
--runThreadN 8 \
--genomeDir /stage/hg19_star_sjdb \
--genomeLoad NoSharedMemory \
--readFilesIn \
EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-unmapped.sorted.fq \ EXAMPLE_PE.rep2_clip.C01.r2.fq.repeat-unmapped.sorted.fq \
--outSAMunmapped Within \
--outFilterMultimapNmax 1 \
--outFilterMultimapScoreRange 1 \
--outFileNamePrefix EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-unmapped.sorted.STAR \
--outSAMattributes All \
--outSAMtype BAM Unsorted \
--outFilterType BySJout \
--outReadsUnmapped Fastx \
--outFilterScoreMin 10 \
--outSAMattrRGline ID:foo \
--outStd Log \
--alignEndsType EndToEnd \
--outBAMcompression 10 \
--outSAMmode Full
# Re-name BAM: rename genome-mapped outputs
mv EXAMPLE_PE.rep2_clip.C01.r1.fq.repeat-unmapped.sorted.STARAligned.out.bam EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mapped.bam
# Name sort BAM: sort output from STAR by name to ensure read pairs are adjacent. samtools sort -n -o EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.bam
EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mapped.bam
## 5. 去除PCR產(chǎn)生的重復(fù)：SE測(cè)序使用umi_tools构拳，常規(guī)工具可以使用barcodecollapsepe.py
# Barcode_collapse_pe (PE): takes output from STAR genome mapping.
barcodecollapsepe.py \
-o EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.rmDup.bam \
-m EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.rmDup.metrics \ -b EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.bam
# Position sort BAM: Takes output from barcode collapse PE (or from SE namesort bam). Sorts resulting bam file for use downstream.
samtools sort -o EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.bam EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDup.bam
# Barcode_collapse_se (SE): takes output from STAR genome mapping. Use umi_tools dedup to identify the extracted random-mer from the previous step and perform PCR duplicate removal.
umi_tools dedup \
--random-seed 1 \
---I EXAMPLE_SE.rep1_clip.umi.r1.fq.genome-mappedSoSo.bam \
--method unique \
--output-stats EXAMPLE_SE.rep1_clip.umi.r1.fq.genome-mappedSoSo.txt \ -S EXAMPLE_SE.rep1_clip.umi.r1.fq.genome-mappedSoSo.rmDup.bam
#Samtools index: Takes output from sortSam, makes bam index for use downstream.
samtools index EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.rmDupSo.bam
## 6. (paired-end only) Merges multiple inline barcodes and filters R1 (uses only R2 for peak calling)
# Samtools merge (PE only): Takes inputs from multiple final bam files. Merges the two technical replicates for further downstream analysis.
samtools merge EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.bam EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.bam EXAMPLE_PE.rep2_clip.D8f.r1.fq.genome- mappedSo.rmDupSo.bam
# Samtools index: Takes output from sortSam, makes bam index for use downstream. samtools index EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.bam
## 7.Calls enriched peak regions (peak clusters) with CLIPPER
# Samtools view (PE only): Takes output from samtools merge. Only outputs the second read in each pair for use with a single stranded peak caller. This is the final bam file to perform analysis on.
# -f 128: 提取read pair的read 2
samtools view -f 128 -b -o EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.bam EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.bam
# Make normalized read density bigwig files. Use --direction f for SE clip as reads are not reversed.
makebigwigfiles \
--bw_pos EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.norm.pos.bw \
--bw_neg EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.norm.neg.bw \
--bam EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.rmDupSo.merged.r2.bam \ --genome hg19.chrom.sizes \
--direction r
# Clipper: Takes results from samtools view. Calls peaks on those files.
clipper \
--species hg19 \
--bam EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.rmDupSo.merged.r2.bam \ --save-pickle \
--outfile EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.peakClusters.bed
## 8. Uses size-matched input sample to normalize and calculate fold-change enrichment within enriched peak regions with custom perl scripts (overlap_peakfi_with_bam_PE.pl, peakscompress.pl)
# Input normalization: Compares the number of reads within the IP sample to the number of reads within the size-matched INPUT sample across Clipper-called peak clusters. This step is performed both within this pipeline as well as within the merge_peaks pipeline using the same perl scripts.
samtools view -cF 4 EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.bam > ip_mapped_readnum.txt samtools view -cF 4 EXAMPLE_PE.rep2_input.NIL.r1.fq.genome- mappedSo.rmDupSo.r2.bam > input_mapped_readnum.txt
overlap_peakfi_with_bam_PE.pl \ EXAMPLE_PE.rep2_clip.C01.r1.fq.genome-mappedSo.rmDupSo.merged.r2.bam \ EXAMPLE_PE.rep2_input.NIL.r1.fq.genome-mappedSo.rmDupSo.r2.bam \ EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.peakClusters.bed \
ip_mapped_readnum.txt \
input_mapped_readnum.txt \
EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.peakClusters.normed.bed
perl compress_l2foldenrpeakfi_for_replicate_overlapping_bedformat.pl \ EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.peakClusters.normed.bed \ EXAMPLE_PE.rep2_clip.C01.r1.fq.genome- mappedSo.rmDupSo.merged.r2.peakClusters.normed.compressed.bed
