师妹的数据线粒体基因组含量30%,FRiP值4%,这个过于低了(标准是0.2)
一方面,用YQY数据做了一下检测流程没问题,他的数据质量比师妹的还要差
另一方面,下载颉伟老师的数据,看一下FRiP值是不是胚胎的就是比较低
1、数据解压
/storage2/anlei/wangwenjing/software/sratoolkit.3.0.2-ubuntu64/bin/fastq-dump --split-3 --gzip SRR2927022.sralite.1
2、TrimGalore过滤
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o /storage2/anlei/MLZ/test/xiewei/clean/ SRR2927022.sralite.1_1.fastq.gz SRR2927022.sralite.1_2.fastq.gz
3、Bowtie2并转格式
bowtie2 -p 5 --very-sensitive -X 2000 -x /storage2/anlei/reference/index/bowtie2/mm10/mm10 -1 SRR2927022.sralite.1_1_val_1.fq.gz -2 SRR2927022.sralite.1_2_val_2.fq.gz -S 8c.sam|samtools sort -@ 5 -O bam -o 8c.bam -
4、去除PCR重复、去除低质量reads
我试了piscard简单一些,步骤比我用samtools少
module load picard/
java -jar $PICARD_PATH MarkDuplicates I=8c.bam O=8c.rmdup.bam M=8c.rmdup.txt REMOVE_DUPLICATES=true
samtools view -h -f 2 -q 10 8c.rmdup.bam | grep -v chrM | samtools sort -O bam -@ 15 -o 8c.last.bam
计算线粒体DNA含量
mtReads=$(samtools idxstats 8c.rmdup.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats 8c.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
老师的是6.92%
!mtDNA
xiewei老师做测序的时候使用了cas9技术消除mtDNA比我们的低很正常
5、Macs2做callpeak
bedtools bamtobed -i 8c.last.bam > 8c.bed
macs2 callpeak -t 8c.bed -g mm --nomodel --shift -100 --extsize 200 -n "8c" --outdir ../peaks/
6、计算FRiP
bed文件是callpeak那步的
Reads=$(bedtools intersect -a 8C.bed -b 8c_peaks.narrowPeak |wc -l|awk '{print $1}')
totalReads=$(wc -l 8c.bed|awk '{print $1}')
echo $Reads $totalReads
echo '==> FRiP value:' $(bc <<< "scale=2;100*$Reads/$totalReads")'%'
老师的是9.53%
!FRiP
这么看来,xiewei老师的数据,即早期胚胎的FRiP值很小,达不到0.2。建议师妹去加测了,做到xiewei老师数据的水平,后续才OK