2023-04-06 xiewei老师8cATAC数据探索

师妹的数据线粒体基因组含量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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容