使用RRBS或WGBS测序数据进行lncRNA基因的甲基化分析

生信小白打算记录一下平时的踩坑过程。

前段时间老师提出一个想法：能不能比较一下早期胚胎发育各个阶段lncRNA和mRNA的甲基化程度差异。我想了一下，如果走完bismark和methylkit整个流程再进行注释，中间需要经历差异分析，这样的话就要选取其中某个阶段作为标准，计算其他阶段相对于此阶段的甲基化程度，其结果不符合预期。因此打算直接从上游分析过程入手，本文将记录整个这一分析流程。

整体思路为：从ensembl网站上下载hg19的lncRNA和cDNA参考基因组，经bismark转化后直接用于测序数据的比对，从比对结果中可以看到总胞嘧啶和甲基化的胞嘧啶的个数，从而得出该样本的甲基化程度，最终剔除掉明显偏离的样本后，取平均值即为该阶段的平均甲基化程度。

首先是前期准备

从GEO数据库下载原始测序数据

prefetch --max-size 50GB --option-file SRR_Acc_List -O ./

将sra文件转换为fastq格式

ls *sra |while read id; do (fastq-dump --split-3 --gzip -O ../lnc/ $id ); done

质控去接头

/home/ngs/anaconda3/bin/trim_galore -q 20 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o ./ SRR${i} _1.fastq.gz SRR${i} _2.fastq.gz #实际运行时此处与后续比对过程为一个小循环#

使用bismark为lncRNA和mRNA参考基因组构建索引

/home/ngs/bin/Bismark-0.22.3/bismark_genome_preparation --bowtie2 --verbose /home/ngs/Documents/library/lncRNA/Homo_sapiens.GRCh38.ncrna.fa
/home/ngs/bin/Bismark-0.22.3/bismark_genome_preparation --bowtie2 --verbose /home/ngs/Documents/library/mRNA/Homo_sapiens.GRCh38.cdna.all.fa

以上两个参考基因组下载于ftp://ftp.ensembl.org/pub/release-101/fasta/homo_sapiens/

然后就是比对过程了，为方便自动化运行，前面的质控和此处一起写了个小循环

for (i=6010333;i<=6010342;i++); do
/home/ngs/bin/Bismark-0.22.3/bismark --bowtie2 --quiet --genome /home/ngs/Documents/library/lncRNA/ --non_directional --fastq -1 SRR${i}_1_val_1.fq,gz -2 SRR${i}_2_val_2.fq,gz -p 10 -o ./
rm SRR${i}_1.fastq.gz SRR${i}_2.fastq.gz SRR${i}_1_val_1.fq,gz SRR${i}_2_val_2.fq,gz #此处是为了腾出些存储空间
done

bismark比对的时间非常久，10线程跑的话，一个样跑一次全部流程大约要4-5小时。等全部样本跑完再写一下根据bismark的比对结果文件计算不同RNA的甲基化程度。