Step1：原始數(shù)據(jù)的準(zhǔn)備
從EMBL數(shù)據(jù)庫(kù)根據(jù)SRR號(hào)下載SRAxxxxxxx.fastq.gz文件（有的是單端測(cè)序，有的是雙末端測(cè)序）
EMBL數(shù)據(jù)庫(kù)網(wǎng)站連接：https://www.ebi.ac.uk/ebisearch/error.ebi
數(shù)據(jù)展示：雙末端：SRR7763560_1.fastq.gz??SRR7763560_2.fastq.gz；
? ? ? ? ? ? ? ? ? ? 單末端：SRR9029562.fastq.gz
數(shù)據(jù)下載linux命令：1.wget? 數(shù)據(jù)的下載連接
? ? ? ? ? ? ? ? 快速下載：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：比對(duì)，利用hisat2講原始數(shù)據(jù)比對(duì)到參考基因組?
1.hisat2下載安裝
2.構(gòu)建hisat2的索引
需要先下載hg38參考基因組，下載hg38.fa.gz文件
文件下載地址：http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz，下載解壓
索引構(gòu)建命令：hisat2-build -p 8 hg38.fa hg38?
大概耗時(shí)半個(gè)小時(shí)，需耐心等待
3.hisat2 比對(duì)命令,生成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文件轉(zhuǎn)換為bam文件
命令：samtools view -bS SRAxxxxxx.sam > SRAxxxxxx.bam
2.去除低質(zhì)量的比對(duì)（MAPQ<30）
命令：samtools view -b -q 30 SRAxxxxxx.bam > SRAxxxxxx.q30.bam
3.利用samtools sort? bam文件得到sort.bam文件，用于接下來(lái)的處理
命令：samtools sort SRAxxxxxx.q30.bam -o?SRAxxxxxx.sort.bam
4.去除duplication
命令：samtools rmdup -s SRAxxxxxx.sort.bam SRAxxxxxx.rmdup.bam
Step4：利用RADAR方法識(shí)別差異甲基化區(qū)域
詳細(xì)參考：https://scottzijiezhang.github.io/RADARmanual/workflow.html#1quick_start
1.將獲得SRAxxxxxx.sort.bam重新命名：
如果是intput數(shù)據(jù)，則重命名為ctrlX.intput.bam【control組】/cX.intput.bam【case組】，如果是IP數(shù)據(jù)則重名為：ctrlX.m6A.bam【control組】/cX.intput.bam【case組】，注意。這里input.bam和m6A.bam中的X是一致的，因?yàn)樗麄兪桥鋵?duì)的。
補(bǔ)充：m6A甲基化數(shù)據(jù)展示：分為control組和case組，每一條數(shù)據(jù)都有配對(duì)的IP和input數(shù)據(jù)，所以重命名時(shí)要很好的區(qū)分IP和input數(shù)據(jù)同時(shí)還要區(qū)分control和case
2.差異分析，需要使用R語(yǔ)言
需要下載hg38對(duì)應(yīng)的注釋文件：
下載地址：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
