教程:b站生信技能树yyds 数据是师妹测完发回来的fq.gz -
DF师兄说了要我建立一个一个文件夹层级储存数据的好习惯
ATAC-seq是双端测序，我现在很大的问题是不会用循环得学，学会写脚本，不能每次都手打代码
第一步：明确样本，构建config.raw（不会写循环此步跳过）
##构建config
ls *_1.fastq.gz > 1
ls *_2.fastq.gz > 2
paste 1 1 2 > config.raw
cat config.raw
第二步：TrimGalore过滤
##trim_galore
-q:设置线程
--phred33 :使用ASCII+33质量值作为Phred得分
--length :去除长度小于参数值的reads
-e :允许的最大误差
--stringency :设置与接头重叠的序列
--paired :对于双端测序文件，正反链质控通过才保留
-o :输出目录
fastq_filel:上一步fastq-dump生成的fastq1文件
fastq_file2:上一步fastq-dump生成的fastq1文件
cat config.raw | while read id;
do echo $id
arr=($id)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o $PATH $fq1 $fq2
#质控 fastqc 一般trimGalore之后都解决的差不多了 就不做了嘻嘻
第三步：Bowtie2比对
#比对后sam转bam
bowtie2 -p 5 --very-sensitive -X 2000 -x /storage2/anlei/reference/index/bowtie2/mm10/mm10 -1 test1.fq -2 test2.fq -S test.sam|samtools sort -@ 5 -O bam -o test.bam -
第四步：去除PCR重复、去除低质量reads
samtools sort -@ 50 -n Aging2.bam -o Aging2.sort.bam
samtools fixmate -@ 50 -m Aging2.sort.bam Aging2.sort.fixed.bam
samtools sort -@ 50 Aging2.sort.fixed.bam -o Aging2.position.fixed.bam
samtools markdup -@ 50 -r Aging2.position.fixed.bam Aging2.rmdup.bam
samtools index Aging2.rmdup.bam
samtools flagstat Aging2.rmdup.bam > Aging2.rmdup.stat
samtools view -h -f 2 -q 30 Aging2.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o - > Aging2.last.bam
samtools index Aging2.last.bam
bedtools bamtobed -i Aging2.last.bam > /storage2/anlei/MLZ/ATAC-seq/callpeak/Aging2.bed
#计算线粒体DNA比例
ps 诺唯赞的技术回的
一般在测序时会注意线粒体DNA的污染率情况。比如少量测序后，线粒体 DNA 的污染率如果没有超过 50%，则不需要进行线粒体 DNA的去除，那么推荐数据量测到 50 M reads/样（15G），以保证有足够的数据量用于后续分析；若线粒体 DNA 的污染率超过 50%，则推荐进行线粒体 DNA 的去除，去除完之后数据量达到 20-30 M reads/样本（6-9G）即可。
## mtReads=$(samtools idxstats Aging2.rmdup.bam | grep 'chrM' | cut -f 3)
## totalReads=$(samtools idxstats Aging2.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
##查看一下过滤情况
##查看两次过滤去除情况
samtools idxstats Aging1.bam | grep -v "_" > 1
samtools idxstats Aging1.rmdup.bam | grep -v "_" > 2
samtools idxstats Aging1.last.bam | grep -v "_" > 3
paste 1 2 3 | cut -f 1,3,7,11
第五步：macs2做call peak
macs2 callpeak -t Aging2.bed -g mm --nomodel --shift -100 --extsize 200 -n "Aging2" --outdir ../peaks/
六、bam转bw
-bamCoverage -p 5 --normalizeUsing CPM -b Aging2.last.bam -o Aging2.bw
七：计算插入片段长度，FRiP，IDR值与deeptools可视化
bamCoverage --normalizeUsing CPM -b Aging1.rmdup.bam -o /storage2/anlei/MLZ/ATAC-seq/bw/Aging1.bw -p 40
可视化
computeMatrix reference-point --referencePoint TSS -p 50 \
-b 10000 -a 10000 \
-R /storage2/anlei/wangwenjing/ChIP-seq/deeptools/bw/mm10_NCBI_refseq.bed \
-S /storage2/anlei/MLZ/ATAC-seq/bw/*.bw \
--skipZeros -o matrix1_TSS.gz \
--outFileSortedRegions regions1_test_genes.bed
1\插入片段
samtools view Aging2.last.bam | cut -f 9 >Aging2.txt
用R画图
A2=read.table('Aging2.txt')
dim(a)
hist(abs(as.numeric(a[,1])), breaks=100)
2、FRiP值的计算
bed文件是callpeak那步的
Reads=$(bedtools intersect -a test_piscard.rmdup.rmchrm.bed -b test_piscard.rmdup.rmchr_peaks.narrowPeak |wc -l|awk '{print $1}')
totalReads=$(wc -l test_piscard.rmdup.rmchrm.bed|awk '{print $1}')
echo $Reads $totalReads
echo '==> FRiP value:' $(bc <<< "scale=2;100*$Reads/$totalReads")'%'
3 、IDR计算
4、相关性计算
multiBamSummary bins --bamfiles /storage2/anlei/MLZ/ATAC-seq/bowtie/cor/*.rmdup.bam --minMappingQuality 30 -out readCounts.npz --outRawCounts readCounts.tab
产生的后缀为npz的文件，通过plotCorrelation命令可以计算相关性，该命令支持spearman和pearson两种相关性分析，pearson相关系数建立在数据符合正态分布的基础上，而spearman相关系数会根据数据的排序即秩进行分析，所以会数据分布没有任何要求，但是对应的敏感性会低一点。
plotCorrelation -in readCounts.npz --corMethod pearson --skipZeros --plotTitle "Pearson Correlation of Average Scores Per Sample" --whatToPlot scatterplot -o scatterplot_PearsonCorr_bigwigScores.pdf --outFileCorMatrix PearsonCorr_bigwigScores.tab
plotCorrelation -in readCounts.npz --corMethod spearman --skipZeros --plotTitle "Sperman Correlation of Read Counts" --whatToPlot heatmap --colorMap RdYlBu --plotNumbers -o heatmap_SpearmanCorr.pdf --outFileCorMatrix SpearmanCorr_readCounts.tab
plotCorrelation -in readCounts.npz --corMethod pearson --skipZeros --plotTitle "pearson Correlation of Read Counts" --whatToPlot heatmap --colorMap RdYlBu --plotNumbers -o heatmap_pearsonCorr.pdf --outFileCorMatrix pearsonCorr_readCounts.tab --removeOutliers
网友评论