fastq测序文件去除rRNA、mapping
input(双端测序)
for i in UV_12h UV_2h WT_2h WT_12h ; do
bowtie2 -x ~/RNASEQ/index/rRNA_index/hg38_rRNA --un-conc-gz ${i}_input_rmrRNA.fastq.gz -1 ${i}_in.read1_Clean.fastq.gz -2 ${i}_in.read2_Clean.fastq.gz -p 8 -S ${i}_rRNA.sam;
rm ${i}_rRNA.sam;
tophat2 -p 8 -G ~/RNASEQ/index/iGenome/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf -o ${i}_out/ ~/RNASEQ/index/iGenome/Homo_sapiens/Ensembl/GRCh37/Sequence/Bowtie2Index/genome ${i}_input_rmrRNA.fastq.1.gz ${i}_input_rmrRNA.fastq.2.gz >${i}tophat2.log;
done
IP (单端测序)
for i in UV_12h UV_2h WT_2h WT_12h ; do
bowtie2 -x ~/RNASEQ/index/rRNA_index/human_rRNA --un-gz ${i}_IP_rmrRNA.fastq.gz -U ${i}_IP.read1_Clean.fastq.gz -p 8 -S ${i}_rRNA.sam; rm ${i}_rRNA.sam;
tophat2 -p 8 -G ~/RNASEQ/index/iGenome/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf -o ${i}_out/ ~/RNASEQ/index/iGenome/Homo_sapiens/Ensembl/GRCh37/Sequence/Bowtie2Index/genome ${i}_ip_rmrRNA.fastq.gz >${i}tophat2.log;
done
MeTPeak calling peak(R)
# using the data included in the package
library(MeTPeak)
# in the real case, change the gtf to what you need
gtf <- "~/RNASEQ/index/iGenome/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf"
ip_uv_2h <- "./ip/UV_12h_out/accepted_hits.bam"
ip_uv_12h <- "./ip/UV_12h_out/accepted_hits.bam"
ip_wt_2h <- "./ip/WT_2h_out/accepted_hits.bam"
ip_wt_12h <- "./ip/WT_12h_out/accepted_hits.bam"
input_uv_2h <- "./input/UV_2h_out/accepted_hits.bam"
input_uv_12h <- "./input/UV_12h_out/accepted_hits.bam"
input_wt_2h <- "./input/WT_2h_out/accepted_hits.bam"
input_wt_12h <- "./input/WT_12h_out/accepted_hits.bam"
metpeak(GENE_ANNO_GTF=gtf,IP_BAM = ip_uv_2h,INPUT_BAM = input_uv_2h,EXPERIMENT_NAME="UV_2H")
metpeak(GENE_ANNO_GTF=gtf,IP_BAM = ip_uv_12h,INPUT_BAM = input_uv_12h,EXPERIMENT_NAME="UV_12H")
metpeak(GENE_ANNO_GTF=gtf,IP_BAM = ip_wt_2h,INPUT_BAM = input_wt_2h,EXPERIMENT_NAME="WT_2H")
metpeak(GENE_ANNO_GTF=gtf,IP_BAM = ip_wt_12h,INPUT_BAM = input_wt_12h,EXPERIMENT_NAME="WT_12H")
网友评论