前面我们分享了 跟着Nature Medicine学MeDIP-seq数据分析,数据和代码都是公开,这个2G的压缩包文件,足以学习3个月,写60篇教程。
这里面的MeDIP-seq指的是DNA,那么MeRIP-seq其实就是RNA水平的又叫做m6a测序,恰好看到了咱们的表观微信交流群我们的生信技能树优秀转录组讲师在分享全套MeRIP-seq文章图表复现代码,我借花献佛整理一下分享给大家:
这个代码关联到了两个 文章,首先是Cell Rep. 2021 Mar ,标题是:《Post-transcriptional regulation of antiviral gene expression by N6-methyladenosine》,然后是文章《Limits in the detection of m6A changes using MeRIP/m6A-seq》
数据在:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE155413 , 是12个样品:
GSM4704026 Ifn input rep1
GSM4704027 Ifn input rep2
GSM4704028 Ifn input rep3
GSM4704029 Ifn meRIP rep1
GSM4704030 Ifn meRIP rep2
GSM4704031 Ifn meRIP rep3
GSM4704032 Mock input rep1
GSM4704033 Mock input rep2
GSM4704034 Mock input rep3
GSM4704035 Mock meRIP rep1
GSM4704036 Mock meRIP rep2
GSM4704037 Mock meRIP rep3
这个实验设计是超级严谨了,其实就Ifn处理前后的两个组,但是每个组都有3个重复,而且是非常严格的MeRIP和 input方便找peaks。
最重要的是,提供了全套MeRIP-seq 图表复现代码在GitHub:https://github.com/al-mcintyre/merip_reanalysis_scripts 这个也是接近2G的压缩包!
首先是文章的重要的图表绘制的R代码 :
2.7K Jan 9 2020 plot_fig1b.R
7.6K Jan 9 2020 plot_fig2.R
20K Jan 9 2020 plot_fig3.R
4.0K Jan 9 2020 plot_fig3_helper.R
11K Jan 9 2020 plot_fig4.R
4.7K Jan 9 2020 plot_fig4_mergedpeaks.R
3.6K Jan 9 2020 plot_fig5.R
6.4K Jan 9 2020 plot_sfig1a.R
4.4K Jan 9 2020 plot_sfig1g.R
1.0K Jan 9 2020 plot_sfig2.R
368B Jan 9 2020 rhelper.R
9.9K Jan 9 2020 run_metdiff_peakcaller.R
然后是具体数据分析的Linux上游的shell脚本:
830B Jan 9 2020 merged_peak_analysis.sh
1.5K Jan 9 2020 plot_fig1a.sh
1.3K Jan 9 2020 plot_fig1bc.sh
3.8K Jan 9 2020 plot_fig2.sh
926B Jan 9 2020 plot_fig4prep.sh
931B Jan 9 2020 plot_sfig2.sh
7.9K Jan 9 2020 run_full_analysis.sh
1.3K Jan 9 2020 run_kallisto.sh
921B Jan 9 2020 run_macs2.2.sh
3.1K Jan 9 2020 run_star.1.sh
研究者自己擅长python,所以也有一些配套的python小脚本:
4.4K Jan 9 2020 get_most_highly_expressed.py
3.4K Jan 9 2020 plot_fig1a_counts_vs_expression.py
799B Jan 9 2020 plot_fig1a_summary.py
1.4K Jan 9 2020 plot_fig1c_replicate_overlap.py
图的颜值不错:
很早以前我就在《生信技能树》发布过教程:新的ngs流程该如何学习(以CUT&Tag 数据处理为例子),提到了我自己是不太可能去把所有的ngs流程全部录制视频的,只能说是更好的传达学习方法给到大家。其实如果你看过我表观组学,比如《ChIP-seq数据分析》 和 《ATAC-seq数据分析》 就会发现其实这个m6A数据处理大同小异的,当然了,肯定是会有一些细微差异是需要注意的。具体数据分析的Linux上游的shell脚本看研究者的GitHub :
7.9K Jan 9 2020 run_full_analysis.sh
1.3K Jan 9 2020 run_kallisto.sh
921B Jan 9 2020 run_macs2.2.sh
3.1K Jan 9 2020 run_star.1.sh
可以看到首先是 trimmomatic 这个java软件进行fastq文件的质量控制:
java -jar trimmomatic-0.36.jar PE -phred33 $PRETRIM1 $PRETRIM2 $FASTQ1 fastqs/${f}_R1_unpaired.fq.gz $FASTQ2 fastqs/${f}_R2_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:15
然后是STAR这个软件进行比对:
STAR --genomeDir $STAR_INDEX --runThreadN $NUM_THREADS --readFilesIn fastqs/${f}_R1.fastq.gz fastqs/${f}_R2.fastq.gz --outFilterMultimapNmax 1 --readFilesCommand zcat --outFileNamePrefix ${f}.star. --genomeLoad LoadAndKeep --outSAMstrandField intronMotif
还有kallisto这个软件进行定量
kallisto quant -t $NUM_THREADS -i $IDX -o kallisto_results/$SAMPLE.kallisto fastqs/${COND}_input_${REP}_R1.trim.fastq.gz fastqs/${COND}_input_${REP}_R2.trim.fastq.gz
到目前为止,其实都是跟转录组流程类似,但是因为有 macs2软件进行 callpeak 代码是:
macs2 callpeak -t alignments/${COND}_IP_${REP}.star.sorted.bam -c alignments/${COND}_input_${REP}.star.sorted.bam --nomodel --extsize $FRAGLEN -g 100e6 -n macs2_results/${COND}_${REP}.macs2 -f BAM --verbose 3
这一切代码能成功运行的前提是Linux基础,再怎么强调生物信息学数据分析学习过程的计算机基础知识的打磨都不为过,我把它粗略的分成基于R语言的统计可视化,以及基于Linux的NGS数据处理:
把R的知识点路线图搞定,如下:
Linux的6个阶段也跨越过去 ,一般来说,每个阶段都需要至少一天以上的学习: