五、Sample QC （发错货这事，有时挺常见的）

除了 早期测序质控，比对后质控，还有一个质控——实验质控（带生物学意义的质控） deeptools真的很deep

1、preparation（两种模式——全基因组|特定区域）

work_dir=/home4/sjshen/project/m6a/GSE29714_Cell
hisat2_bam_dir=$work_dir/1.1_hisat2/mouse/position_sorted_bam

echo $hisat2_bam_dir
cd $hisat2_bam_dir

(nohup multiBamSummary bins  -p 60 --bamfiles $hisat2_bam_dir/*bam -out ./deeptools/results_whole_region.npz ) &
(nohup multiBamSummary BED-file --BED ./mm10.RefSeq.bed  -p 60 --bamfiles $hisat2_bam_dir/*bam -out ./deeptools_qc/results_gene_region.npz ) &

Tips1：multiBamSummary是deeptools套件的一个函数，主要用于生成下游质控的初始文件 Tips2：默认划bins是10kb一个bin Tips3：mm10.RefSeq.bed 来自 UCSC TABLE broswer

2、Correlation

plotCorrelation -in results_gene_region.npz --corMethod spearman --skipZeros --plotTitle "Spearman Correlation of Read Counts gene region" \
--labels SYSY-MeRIP-rep1 input-rep1 NEB-MeRIP-rep2 SYSY-MeRIP-rep2 input-rep2 --colorMap PiYG \
--whatToPlot heatmap  --plotNumbers -o heatmap_SpearmanCorr_readCounts_gene_region.png --outFileCorMatrix SpearmanCorr_readCounts_gene_region.tab &

plotCorrelation -in results_gene_region.npz --corMethod spearman --skipZeros --plotTitle "Spearman Correlation of Read Counts gene region" \
--labels SYSY-MeRIP-rep1 input-rep1 NEB-MeRIP-rep2 SYSY-MeRIP-rep2 input-rep2 --log1p \
--whatToPlot scatterplot -o heatmap_SpearmanCorr_readCounts_gene_region_2.png &

Tips1：spearman 和 pearson 是最常用的方法，pearson更容易受异常值影响 Tips2：用于判断异常样本 Tips3：文件名太长可以用 –labels 修改 Tips4：可修改颜色

heatmap

scatterplot

3、PCA

cd  $hisat2_bam_dir/deeptools_qc

plotPCA --transpose --markers 'o' '>' 'o' 'o' '>' -in results_gene_region.npz -o PCA_readCounts_gene_region.png -T "PCA of read counts gene region" --outFileNameData PCA_readCounts_data_gene_region.tab &

plotPCA --transpose --markers 'o' '>' 'o' 'o' '>' -in results_whole_region.npz -o PCA_readCounts_whole_region.png -T "PCA of read counts whole region" --outFileNameData PCA_readCounts_data_whole_region.tab &

Tips1：PCA是一种降维方法，他将样本的诸多差异（基因|bins|…）浓缩到二维，其中PC1和PC2是区别最大的前两个维度 Tips2：默认取前1000个影响因素 Tips3：用于判断异常样本 Tips4：文件名太长可以用 –labels 修改

gene_region范围计算

whole_region范围计算

4、fragment length

work_dir=/home4/sjshen/project/m6a/GSE29714_Cell
hisat2_bam_dir=$work_dir/1.1_hisat2/mouse/position_sorted_bam

for i in $hisat2_bam_dir/*.bam
do
echo $i
(nohup bamPEFragmentSize -p 20 -hist ${i##*/}fragmentSize.png -T "Fragment size of SE MeIP-seq data" --maxFragmentLength 200 --bamfiles  $i  --samplesLabel  ${i##*/}) > ./${i##*/}fragmentSize.log &
done

Tips：主要用于判断片段超声情况，一般m6A大约在100-200nt

5、other

plotFingerprint
plotCoverage
#这两个函数主要检查抗体富集和reads覆盖度，仅供参考