1.下載數(shù)據(jù)：ascp/prefetch。我們的服務(wù)器不要用ascp下，不要用ascp下，不要用ascp下?。?！下不下來(lái)的。
命令：prefetch.2.9.3 -v SRR7366879。下載的數(shù)據(jù)自動(dòng)存在~/ncbi/public/sra文件夾下，這個(gè)文件夾塞滿(mǎn)了會(huì)用不了tab鍵。移動(dòng)到自己的文件夾。
用fastq-dump.2.9.6 --gzip --split-3 SRR7366882.sra -O Dir/RIPseq/SRR7366882將sra數(shù)據(jù)轉(zhuǎn)換為fastq格式。
因?yàn)橹挥?個(gè)文件夾，所以是一個(gè)一個(gè)下的，如果數(shù)據(jù)比較多，可以寫(xiě)個(gè)循環(huán)
2.質(zhì)控：fastqc
命令：fastqc -o dir/RIPseq/SRR7366882/fastqc -t 12 dir/RIPseq/SRR7366882/SRR7366882_1.fastq
考慮去除前4個(gè)堿基。另一個(gè)是Kmer Content也未過(guò)檢測(cè),暫時(shí)不管。采用trimmomatic質(zhì)控之后質(zhì)量合格。trimmomatic參數(shù)：java -jar /home/softwares/Trimmomatic-0.33/trimmomatic-0.33.jar PE -threads 20 dir/SRR7366882_1.fastq dir/SRR7366882_2.fastq dir/SRR7366882_1P.fq dir/SRR7366882_1U.fq dir/SRR7366882_2P.fq dir/SRR7366882_2U.fq ILLUMINACLIP:/home/softwares/Trimmomatic-0.33/adapters/TruSeq3-PE.fa:2:30:12:1:true HEADCROP:4 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
3.比對(duì)：star
第一步：star產(chǎn)生對(duì)基因組產(chǎn)生index:
star比對(duì)得到bam文件，對(duì)得到的bam文件進(jìn)行質(zhì)控：
首先得到Uniq的reads:
samtools view -H SRR7366882.bam>SRR7366882_header
samtools view SRR7366882_mapped_reads.bam|grep -w "NH:i:1"|cat SRR7366882_header -|samtools view -b>SRR7366882_mapped_reads_uniq.bam &
再得到去除Pcr重復(fù)之后的reads:
samtools rmdup ${dir}"_mapped_reads_uniq.bam" ${dir}"_mapped_reads_rmdup.bam"
4.peak calling:MACS2
/mnt/data/slt/lq/RIPseq/alignment/SRR7366880/SRR7366880_mapped_reads.bam?-g hs --nomodel --extsize 200 -n tumor1 -f BAM--outdir /mnt/data/slt/lq/RIPseq/peak_calling/tumor1/tumor1 &
文章說(shuō)基于RNA片段的長(zhǎng)度設(shè)置--shift 200，可是我覺(jué)得這有問(wèn)題，因?yàn)榘凑誱acs方法文章的說(shuō)法，shift應(yīng)該是絕對(duì)偏移量。macs2本來(lái)是為了call轉(zhuǎn)錄因子結(jié)合的峰，由于實(shí)際上測(cè)不到轉(zhuǎn)錄因子的結(jié)合區(qū)域，所以需要把seq數(shù)據(jù)偏移一定距離以更好的得到轉(zhuǎn)錄因子結(jié)合的峰?，F(xiàn)在我們把它照搬到m6Aseq上來(lái)，由于m6A本身只是一個(gè)位點(diǎn)，所以我認(rèn)為是不需要這個(gè)偏移量，直接設(shè)置extsize 200就可以的。所以按照我的理解這個(gè)參數(shù) 按照文章中的參數(shù)命令行是macs2?callpeak-t DIR/SRR7366879_mapped_reads.bam -c DIR/SRR7366880_mapped_reads.bam -g hs --nomodel?--extsize?200 --shift 100 -n tumor1 -f BAM--outdir DIR/tumor1。不過(guò)一般不用指定--extsize 200，--shift 100采用自動(dòng)建模，除非明確知道RNA fragment長(zhǎng)度。
5.Differential m6A level analysis
首先將2個(gè)腫瘤中narrow peak有重疊的區(qū)域合并，沒(méi)有重疊的也保留（bedtools merge），合起來(lái)作為一個(gè)reference m6A peaks list，然后用bedtools 的multicov工具去計(jì)算input組和tumor組在這個(gè)區(qū)域的count數(shù)，fisher檢驗(yàn)差異是否顯著。分別計(jì)算兩個(gè)tumor相對(duì)于input的count ratio，得到log2[RTumor1/RTumor2],>1表示在tumor1中m6A level水平較高，數(shù)據(jù)如下圖所示。
可根據(jù)p_values篩選出IP組特異的峰，根據(jù)log2[RTumor1/RTumor2]進(jìn)一步選出特定腫瘤較高的峰。加上uniq的數(shù)據(jù)可以簡(jiǎn)單統(tǒng)計(jì)：
問(wèn)題：文章里說(shuō)Differential?m6A level analysis是用的count數(shù)據(jù)，我也是用的count數(shù)，都是整數(shù)。但是附加文件中的raw data包含小數(shù)：
另外作者畫(huà)了這個(gè)圖（s6B）,原始數(shù)據(jù)(文章S6B_rawdata),不太懂其中l(wèi)og2(IP_FC)怎么來(lái)的
