十、差异分析（五指兄弟不同长，天地万物终有别）

我们实验总有 对照组/处理组 或者 时间线，所以多样本的分析也很重要

1、peak related gene 交并补-差异 （Venn来了，peak定位好时推荐）

《=== Venn示例 ===》 ###请打开 Rcode - part10_1 《=== Venn神器 ===》

2、difference based on value per region (peak定位不好时，可以以region来做差异)

# Step1 - preparation -  region可以来源于现有peak，也可以全基因组计算(bedtools makewindows)

cd  /home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam/peak_operation
work_dir=/home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam/peak_operation

# intersect
bedtools intersect -a $work_dir/peak_bed_home/SYSY_rep1_peaks.narrowPeak  -b $work_dir/peak_bed_home/SYSY_rep2_peaks.narrowPeak > $work_dir/intersect/SYSY_common.bed
bedtools intersect -a $work_dir/peak_bed_home/NEB_rep2_peaks.narrowPeak  -b $work_dir/peak_bed_home/SYSY_rep2_peaks.narrowPeak  > $work_dir/intersect/NEB_SYSY_common.bed

# merge (cat/sort/merge)
cat $work_dir/peak_bed_home/SYSY_rep1_peaks.narrowPeak $work_dir/peak_bed_home/SYSY_rep2_peaks.narrowPeak $work_dir/peak_bed_home/NEB_rep2_peaks.narrowPeak > $work_dir/merge/all_common.bed
bedSort all_common.bed all_common.sorted.bed
bedtools merge -i all_common.sorted.bed > all_common.merge.bed

#--------------------------------------------------------------------

# Step2 - calculate peak RPKM 

cd  /home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam/peak_operation
work_dir=/home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam/peak_operation

#  A - count
#bed_dir=/home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam/peak_operation/merge
#bed_dir=/home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam/peak_operation/intersect

bam_dir=/home4/sjshen/project/m6a/GSE29714_Cell/1.1_hisat2/mouse/position_sorted_bam

# cd $bam_dir
# for i in *.bam
# do
# echo $i
# samtools flagstat -@ 10 $i > ./${i/.bam/.flagstat} &
# done

for i in $bed_dir/*.bed
do
t1=${i##*/}
bed_name=${t1%%.bed}
echo $bed_name
    for j in $bam_dir/*.bam
    do
    t2=${j##*/}
    bam_name=${t2%%.hisat*}
    echo $bam_name
    (nohup multiBamSummary BED-file --BED $i -p 20 --centerReads --bamfiles $j --outRawCounts ./${bed_name}.${bam_name}.readCounts.bed -out ./${bed_name}.${bam_name}.readCounts.npz )&
    done
done

#  B - RPKM
for i in *.readCounts.bed
do
t1=${i%%.readCounts.bed}
name=${t1##*.}
echo $name
total_mapped_hit=`awk ' BEGIN{FS="[ \t();%]+";OFS="\t"} FNR==5 {print $1}' $bam_dir/${name}.hisat2.position_sorted.flagstat`
echo $total_mapped_hit
awk -v total_mapped_hit=$total_mapped_hit 'BEGIN{FS="\t";OFS="\t"} NR>1 {print $1,$2,$3,($4*1e6)/(total_mapped_hit*($3-$2))*1000} ' $i > ${t1}.deeptools.rpkm.bed &
done

#total_mapped_hit=`awk ' BEGIN{FS="[ \t();%]+";OFS="\t"} FNR==9 {print $1}' $bam_dir/${name}.flagstat`

# C - sort
for i in *.deeptools.rpkm.bed
do
echo $i
bedSort $i ${i/rpkm/rpkm.sorted} &  ###ucsc tools
done

# D - clean 
for i in *.deeptools.rpkm.sorted.bed
do
echo $i
egrep -v \(Pf3D7_API_v3\|Pf_M76611\) $i > ${i/.deeptools.rpkm.sorted.bed/.rpkm.clean.bed} &
done

#--------------------------------------------------------------------

# Step3 - plot in R
# 详见代码
>

  《=== 示例备份 ===》 ###请打开 Rcode - part10_2

Dependencies

• Softwares
  • samtools
  • BEDtools
  • deeptools