1.下载数据：ascp/prefetch。我们的服务器不要用ascp下，不要用ascp下，不要用ascp下！！！下不下来的。
命令：prefetch.2.9.3 -v SRR7366879。下载的数据自动存在~/ncbi/public/sra文件夹下，这个文件夹塞满了会用不了tab键。移动到自己的文件夹。
用fastq-dump.2.9.6 --gzip --split-3 SRR7366882.sra -O Dir/RIPseq/SRR7366882将sra数据转换为fastq格式。
因为只有4个文件夹，所以是一个一个下的，如果数据比较多，可以写个循环
2.质控：fastqc
命令：fastqc -o dir/RIPseq/SRR7366882/fastqc -t 12 dir/RIPseq/SRR7366882/SRR7366882_1.fastq
考虑去除前4个碱基。另一个是Kmer Content也未过检测,暂时不管。采用trimmomatic质控之后质量合格。trimmomatic参数：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.比对：star
第一步：star产生对基因组产生index:
star比对得到bam文件，对得到的bam文件进行质控：
首先得到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重复之后的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 &
文章说基于RNA片段的长度设置--shift 200，可是我觉得这有问题，因为按照macs方法文章的说法，shift应该是绝对偏移量。macs2本来是为了call转录因子结合的峰，由于实际上测不到转录因子的结合区域，所以需要把seq数据偏移一定距离以更好的得到转录因子结合的峰。现在我们把它照搬到m6Aseq上来，由于m6A本身只是一个位点，所以我认为是不需要这个偏移量，直接设置extsize 200就可以的。所以按照我的理解这个参数 按照文章中的参数命令行是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。不过一般不用指定--extsize 200，--shift 100采用自动建模，除非明确知道RNA fragment长度。
5.Differential m6A level analysis
首先将2个肿瘤中narrow peak有重叠的区域合并，没有重叠的也保留（bedtools merge），合起来作为一个reference m6A peaks list，然后用bedtools 的multicov工具去计算input组和tumor组在这个区域的count数，fisher检验差异是否显著。分别计算两个tumor相对于input的count ratio，得到log2[RTumor1/RTumor2],>1表示在tumor1中m6A level水平较高，数据如下图所示。
可根据p_values筛选出IP组特异的峰，根据log2[RTumor1/RTumor2]进一步选出特定肿瘤较高的峰。加上uniq的数据可以简单统计：
问题：文章里说Differential m6A level analysis是用的count数据，我也是用的count数，都是整数。但是附加文件中的raw data包含小数：
另外作者画了这个图（s6B）,原始数据(文章S6B_rawdata),不太懂其中log2(IP_FC)怎么来的
