Step1：原始数据的准备
从EMBL数据库根据SRR号下载SRAxxxxxxx.fastq.gz文件（有的是单端测序，有的是双末端测序）
EMBL数据库网站连接：https://www.ebi.ac.uk/ebisearch/error.ebi
数据展示：双末端：SRR7763560_1.fastq.gz SRR7763560_2.fastq.gz；
单末端：SRR9029562.fastq.gz
数据下载linux命令：1.wget 数据的下载连接
快速下载：2.ascp-QT-l20m-P33001-i /home/wangi/anaconda2/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR110/000/SRR1104940/SRR1104940.fastq.gz /home/wanghongli/histone/sheet46/GSE53938/SRR1104940.fastq.gz
Step2：比对，利用hisat2讲原始数据比对到参考基因组
1.hisat2下载安装
2.构建hisat2的索引
需要先下载hg38参考基因组，下载hg38.fa.gz文件
文件下载地址：http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz，下载解压
索引构建命令：hisat2-build -p 8 hg38.fa hg38
大概耗时半个小时，需耐心等待
3.hisat2 比对命令,生成sam文件
双末端：
hisat2 -p 16 -t -x /hisat2索引所在的文件路径/hg38 -1 SRR7763560_1.fastq.gz -2 SRR7763560_2.fastq.gz -S SRR7763560.sam
单末端：
hisat2 -p 16 -t -x /hisat2索引所在的文件路径/hg38 -U SRR9029562.fastq.gz -S SRR9029562.sam
Step3：
1.利用samtools将sam文件转换为bam文件
命令：samtools view -bS SRAxxxxxx.sam > SRAxxxxxx.bam
2.去除低质量的比对（MAPQ<30）
命令：samtools view -b -q 30 SRAxxxxxx.bam > SRAxxxxxx.q30.bam
3.利用samtools sort bam文件得到sort.bam文件，用于接下来的处理
命令：samtools sort SRAxxxxxx.q30.bam -o SRAxxxxxx.sort.bam
4.去除duplication
命令：samtools rmdup -s SRAxxxxxx.sort.bam SRAxxxxxx.rmdup.bam
Step4：利用RADAR方法识别差异甲基化区域
详细参考：https://scottzijiezhang.github.io/RADARmanual/workflow.html#1quick_start
1.将获得SRAxxxxxx.sort.bam重新命名：
如果是intput数据，则重命名为ctrlX.intput.bam【control组】/cX.intput.bam【case组】，如果是IP数据则重名为：ctrlX.m6A.bam【control组】/cX.intput.bam【case组】，注意。这里input.bam和m6A.bam中的X是一致的，因为他们是配对的。
补充：m6A甲基化数据展示：分为control组和case组，每一条数据都有配对的IP和input数据，所以重命名时要很好的区分IP和input数据同时还要区分control和case
2.差异分析，需要使用R语言
需要下载hg38对应的注释文件：
下载地址：ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.annotation.gtf.gz
Rcode:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("devtools")
library(devtools)
install_github("RADAR")
library(RADAR)
radar <- countReads(
samplenames = c("ctrl1","ctrl2","ctrl3","c1","c2","c3"),
gtf = "/注释文件所在文件路径/gencode.v32.annotation.gtf",
bamFolder = "bam文件所在的文件路径",
modification = "m6A",
outputDir = "文件输出路径",
threads = 10
)
radar <- normalizeLibrary( radar )
radar <- adjustExprLevel( radar )
variable(radar) <- data.frame( Group = c("Ctl","Ctl","Ctl","Case"，"Case"，"Case"，) )
radar <- filterBins( radar ,minCountsCutOff = 15)
radar <- diffIP_parallel(radar,thread = 10)
radar <- reportResult( radar, cutoff = 0.1, Beta_cutoff = 0.5 )
res <- results( radar )
write.table(res,file='/home/wanghongli/m6adata/GSE87190/MA9.3ITD/0911_diff.xls',sep='\t',quote=F)
samtools view -b -q 30 $i.bam > $i.q30.bam
