写在前面

可能第一眼看到这个标题大家都会很奇怪，因为2G内存的服务器还要分析上游的话，不太可能，确实在中间我把比对的过程跳过了
服务器生信分析平台配置参考https://zhuanlan.zhihu.com/p/74423987
生信是实践性很强的交叉学科，通过实践提高自己思考和解决问题的能力

本文结合linux操作和R，R单独在本地其实也可以做，用esATAC，一步包装到位就是

主要参考文章是生信技能树上的教程：https://www.jianshu.com/p/5bce43a537fd

和：https://mp.weixin.qq.com/s/a4qAcKE1DoukpLVV_ybobA

然而这个包在我电脑上run报错了：

Error: pandoc document conversion failed with error 11
In addition: There were 50 or more warnings (use warnings() to see the first 50)

改天测试

current work plan：

• RNAseq下游分析在R里面比较简单，上下游结合的教程整理；
• TCGA 数据库挖掘，推荐在cellpress上的泛癌研究文献：https://www.cell.com/pb-assets/consortium/pancanceratlas/pancani3/index.html
• GEO-master的airwayRNA
• ChIPpeakAnno, atacseq proj分析报告
• 加拿大生物信息学研讨会资源宝藏：https://bioinformatics.ca/workshops/2018-epigenomic-data-analysis/
  • 对应youtube链接：https://www.youtube.com/channel/UCKbkfKk65PZyRCzUwXOJung
• 单细胞培训：https://www.csc.fi/web/training/-/scrnaseq可以下载code

在win10上用linux子系统做RNAseq

参考的视频：https://www.bilibili.com/video/av65970178?p=4

参考笔记和paper: https://www.jianshu.com/p/6fa6da3bae40

胆子够大可以尝试搭建linux虚拟机

下面是几个points：

1. win10默认挂载到/mnt上，这对于拷贝文件很方便，不用FileZilla传来传去

2. 下载的软件zip、安装的软件路径、编译软件等目录要安排好,养成好习惯，conda配置环境work时也是如此:

   *bin：可执行程序，添加所有软件链接方便加入环境变量*
   
   *biosoft：软件安装目录，解压编译后放这里*
   
   *src：软件源代码，最开始下载的包放在这里*
   

   下载软件失败，发现原来在备份sources.list的时候把sources打成了source

3. 操作.bashrc文件重新个性化配置时（用户的.bashrc里放着用户自己的环境变量,export $PATH可以添加软件路径到里面，推荐用conda解决大部分软件环境问题），vim的操作要熟悉

4. 编译gffcompare源代码的时候，g++编译器下载失败，apt-get做不了，说是有一系列依赖包没下，然后我发现那些依赖包也下不了，想放弃的时候从头检查感觉可能是etc/apt下的sources.list文件有问题，把阿里云的软件下载镜像站点都换成清华的以后，编译gffcompare就成功了

5. HISAT2速度相当快,stringtie重构转录本同时组装转录本和评估其表达水平，gffcompare评估组装效果（检查拼接的转录本有多少匹配到了参考基因组注释文件，gffcompare可以用来compare, merge, annotate and estimate accuracy of one or more GFF files), 最后再Stringtie估算转录本表达量并生成Ballgown能使用的统计表格，后面就是下游的事了

   
   
   stringtie.png

如上图，stringtie是新转录本，可能是lncRNA

偶然间看到腾讯云特价的学生用服务器后果断买了，最便宜的1核2G学生用服务器，尽管远远不够生信分析的要求，体验还是可以的，有条件的同学建议买16核64G服务器，如果想跑很多个样本的数据的话（比如单细胞）。

下载数据折腾3-4天

因为网速/有的站点网络不稳定的问题，下载数据很多时候也很麻烦

国内站点下载数据一直是个问题，后来才看到有人推荐用海外云服务器，按用量计费（可能比较划算），手动下载是不现实的（不支持用wget，用axel多线程断点续传花了两天才下了一个3G的测试数据（HISAT文件的RNAseq pipeline里的）)听说aspera超高速下载sra数据很棒，尝试了一下失败（ncbi可能取消了ascp协议，或者我自己配置的问题，不支持ncbi下载，还有一个重要原因是ncbi的ftp目录组织不太好，经常会变化）但是ncbi不是唯一的出路

最后选择了EBI-ERA数据库下载，同样支持SRR号下载，而且界面科学，点击ftp复制链接即可，而且下载后直接是fasta文件，省去了ncbi下载的sra文件还要fastq-dump的步骤，而且ERA是绝对支持ascp的！速度超快几百M/s，即使数据量很大但也没花多少时间，这个时候考虑到我的硬盘要省着点用（因为后面测试时因存储问题，文件是truncated的，只好删掉两个测试数据）

#注意密钥文件位置，vol1后根据自己要下载数据的SRR号修改url即可
ascp -v -k 1 -T -l 200m -P 33001 -i ~/miniconda3/envs/atac/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1

下载数据在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581

conda配置环境

照搬代码，不解释，注意是py2，安装完一定要检查软件有没有安装成功，否则像我一样后面又删除环境重装。添加channel后看一下channel的优先级，建议把bioconda放最上面，减少search时间。

关于conda的基础操作：https://www.jianshu.com/p/edaa744ea47d

推荐安装软件的课程：https://www.bilibili.com/video/av82983858?p=7

所有的软件安装失败都是环境配置造成的

# https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
# https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/ 
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc 
## 安装好conda后需要设置镜像。
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

conda  create -n atac -y   python=2 bwa
conda info --envs
source activate atac
# 可以用search先进行检索
conda search trim_galore
## 保证所有的软件都是安装在atac这个环境下面
conda install -y sra-tools  
conda install -y trim-galore  samtools bedtools
conda install -y deeptools homer  meme
conda install -y macs2 bowtie bowtie2 
conda install -y  multiqc 
conda install -y  sambamba

(atac) msq 20:54:13 ~
$ ll
total 28
drwxrwxr-x  2 msq msq 4096 Feb 28 14:46 bin
drwxrwxr-x  2 msq msq 4096 Feb 28 14:46 biosoft
drwxrwxr-x 16 msq msq 4096 Feb 28 13:32 miniconda3
drwxrwxr-x  3 msq msq 4096 Feb 29 12:34 project
drwxr-xr-x  3 msq msq 4096 Feb 28 08:44 R
-rw-rw-r--  1 msq msq  100 Feb 29 20:08 removelog.sh
drwxrwxr-x  2 msq msq 4096 Feb 28 14:19 src

上面这个是我的home目录，用来做atac等流程的目录均放在project下。建议在home下放置一个sh文件专门用来清空log文件，否则时间长了特别占存储

质控

fastqc + multiqc两件套，不多解释，可以好好看看输出的报告文件，没安装的自己安装

参考：https://www.baidu.com/link?url=oukLwbl54rtq3gFEAbR05ucV6HU8K6pnb-X1IWLaR6d3P4tSWWH0cxptT61r_pbi&wd=&eqid=f427e0fd00153c79000000065e5e43c9

fastqc官网，含示例报告文件：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

推荐trim_galore和shell编程

trim_galore去接头，这个软件非常智能，自动探测adaptor，比trimmomatic好

参考该软件：https://www.jianshu.com/p/7a3de6b8e503

参考优秀的shell笔记：https://www.jianshu.com/nb/13897796

linux命令查找说明书: https://man.linuxde.net/

cd ~/project/atac/
mkdir -p clean 
source activate atac  
# trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o clean/ raw/2-cell-1_1.fastq.gz raw/2-cell-1_2.fastq.gz
#cat config.raw  |while read id;
#do echo $id
#arr=($id)
#fq2=${arr[2]}
#fq1=${arr[1]}
#sample=${arr[0]}
#nohup  trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o  clean  $fq1   $fq2  & 
#done 
ps -ef |grep trim

当时运行的界面开头是一些运行参数，可以看到AUTO-DETECTING adaptor是Nextera而非illumina，这就是选择trim_galore避免踩坑的原因
trim完以后重新fastqc一下，可以对比：

cd ~/project/atac/qc 
mkdir -p clean 
fastqc -t 5  ../clean/2-cell-*gz -o clean 
mkdir -p raw 
fastqc -t 5  ../raw/2-cell-*gz -o clean 
# https://zh.wikipedia.org/wiki/ASCII
## 还有很多其它工具，比如：
qualimap='/home/jianmingzeng/biosoft/Qualimap/qualimap_v2.2.1/qualimap'
$qualimap bamqc --java-mem-size=20G  -bam $id  -outdir ./ 
#注意multiqc最好在当前都是fastqc report的路径下
multiqc .

比对

比对计算量相当大，这里我的服务器超内存了：processing memory不足以运行bowtie2（忘记截屏）

wget下载索引还是比较慢，最好不要自己构建index

# 索引大小为3.2GB， 不建议自己下载基因组构建，可以直接下载索引文件，代码如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip

因为我自己服务器没法run bowtie（bowtie2对人类全基因组的比对至少消耗4G以上吧），但代码还是可以写的：

#把sort直接一步做了
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq  -2 test2.fq  |samtools sort -@ 5 -O bam -o test.bam -

#计划下载bam文件，直接往下做,index要写全路径（第一次测试的bam文件是在esATAC中的Example.bam，但run到后#面我发现call peak文件是空的，可能是这个文件问题，所以我自己另外下载了bam进行测试)
(atac) msq 21:23:47 ~/project/atac/reference
$ ll
total 20642232
-rw-rw-r-- 1 msq msq 5582057939 Mar  3 17:05 D_trm.bam
-rw-rw-r-- 1 msq msq 4324628903 Mar  3 19:16 D_trm.sort.bam
drwxrwxr-x 3 msq msq       4096 Mar  3 11:06 index #index是我第一次测试时用来放mm10 index的
-rw-rw-r-- 1 msq msq 6520916926 Mar  3 17:42 P_trm.bam
-rw-rw-r-- 1 msq msq 4710012734 Mar  3 21:23 P_trm.sort.bam
#因为文件较大，sort会花些时间

探索sam文件

（或者用samtools view看bam文件），有详细的reads mapping信息，mapping分数

[SAMtools] 常用指令总结

理解并操作bam文件：

很多NGS文件都分成header和record部分，header里面还会包括此文件生成的过程信息,没有sort的bam文件在开头是这样的：

@HD VN:1.0 SO:unsorted

sort后：

@HD VN:1.0 SO:coordinate

header内容不多，也不会太复杂，每一行都用‘@’ 符号开头，里面主要包含了版本信息，序列比对的参考序列信息，如果是标准工具（bwa，bowtie，picard）生成的BAM，一般还会包含生成该份文件的参数信息（如上图），@PG标签开头的那些。这里需要重点提一下的是header中的@RG也就是Read group信息，这是在做后续数据分析时专门用于区分不同样本的重要信息。它的重要性还体现在，如果原来样本的测序深度比较深，一般会按照不同的lane分开比对，最后再合并在一起，那么这个时候你会在这个BAM文件中看到有多个RG，里面记录了不同的lane，甚至测序文库的信息，唯一不变的一定是SM的sample信息，这样合并后才能正确处理。

所以如果一个处理组有好几个重复，还需要把这些bam文件merge一下：samtools merge一个output三个input

写循环除了while read id do还可以用xargs命令：

ls *.bam |xargs -i  samtools index {}

参考链接：https://www.jianshu.com/p/364e640d3c9f

可以对sort后的文件构建index，方便查看的时候直接跳到自己感兴趣的区域查看：

samtools index in.bam  # 生成in.bam的索引文件in.bam.bai
samtools view in.bam chr22            # 跳转到chr22染色体
samtools view in.bam chr22:16050103   # 跳转到chr22:16050103位置
samtools view in.bam chr22:16050103-16050103  # 只查看该位置

另外发现有一些比较奇特的染色体名：（我用的是mm10,mm是20对染色体)

@SQ     SN:chr1_GL456210_random LN:169725
@SQ     SN:chr1_GL456211_random LN:241735
@SQ     SN:chr1_GL456212_random LN:153618
@SQ     SN:chr1_GL456213_random LN:39340
@SQ     SN:chr1_GL456221_random LN:206961

@SQ     SN:chrUn_GL456239       LN:40056
@SQ     SN:chrUn_GL456367       LN:42057
@SQ     SN:chrUn_GL456378       LN:31602
@SQ     SN:chrUn_GL456381       LN:25871
@SQ     SN:chrUn_GL456382       LN:23158

查教程发现：

“chr*_random sequences” 知道来自哪条染色体但不知道具体位置的序列

“chrUn_* sequences” 知道来自人类基因组序列，但不知道与染色体的关系

详细参考：https://www.jianshu.com/p/2cd5e6e310c9

还可以对bam/sam文件做统计：

samtools flagstat D_trm.sort.bam -O D_trm.sort.tsv

输出结果：

##总的reads
62378174 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
#总体reads比对率
49140625 + 0 mapped (78.78% : N/A)
62378174 + 0 paired in sequencing
31189087 + 0 read1
31189087 + 0 read2
#PE reads比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
35840362 + 0 properly paired (57.46% : N/A)
#PE中两条都比对到参考序列上的reads数
47374920 + 0 with itself and mate mapped
#单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数
1765705 + 0 singletons (2.83% : N/A)
#PE中两条分别比对到两条不同的参考序列的reads数
1226554 + 0 with mate mapped to a different chr
#PE中两条分别比对到两条不同的参考序列的reads数（比对质量>=5）
19944 + 0 with mate mapped to a different chr (mapQ>=5)

过滤BAM文件

ATAC-Seq与其他方法不同的一点是需要过滤去除线粒体（如果是植物，还需要过滤叶绿体），因为线粒体DNA是裸露的，也可以被Tn5酶识别切割

参考：https://www.jianshu.com/p/5bddf9c9ca4c

过滤唯一比对的reads

#这句话-h指定打印头文件信息，-t指定线程，-f指定bam文件格式，-F指定过滤条件，[XS]是sam的可选域，描述序列多重比对上时第二个最有比对的得分，只有当序列比对到多个位置时才会出现，因此[XS] == null除去了多重比对的序列，not unmapped 和not duplicate分别用来除去未比对上序列和重复序列
sambamba view -h -t 2 -f bam \
-F "[XS] == null and not unmapped " \
Example.sort.bam > Example_aln.bam

#可以看到过滤以后，总体比对率达到了100%
32526963 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
32526963 + 0 mapped (100.00% : N/A)
32526963 + 0 paired in sequencing
16455894 + 0 read1
16071069 + 0 read2
22927379 + 0 properly paired (70.49% : N/A)
31159503 + 0 with itself and mate mapped
1367460 + 0 singletons (4.20% : N/A)
29462 + 0 with mate mapped to a different chr
7025 + 0 with mate mapped to a different chr (mapQ>=5)

去除PCR重复

这个我没有做，因为上面数据里没有duplicates

java -jar picard-tools-1.119/MarkDuplicates.jar REMOVE_DUPLICATES=true I=H1hesc_Input_Rep1_chr12_aln.bam O=H1hesc_Input_Rep1_chr12_aln.dedup.bam M=H1hesc.duplicates.log

过滤线粒体基因

samtools view -h Example_aln.bam | grep -v 'chrM' | samtools view -bS -o Example.final.bam
#grep -v参数反选
#此时就得到了唯一比对且已经去除过线粒体的比对文件，可以用于接下来的peaks calling

将bam文件转为bed文件

bedtools bamtobed -i Example.final.bam > Example.bed

或者直接这样按照mapq值筛选，-q就是按照mapq筛选

samtools view -bhS -q 30 input.sam > output.bam

最后得到final文件：

$ ll
total 24306500
-rw-rw-r-- 1 msq msq 2033637104 Mar  4 10:46 D.final.bam
-rw-rw-r-- 1 msq msq 5582057939 Mar  3 17:05 D_trm.bam
-rw-rw-r-- 1 msq msq 4324628903 Mar  3 19:16 D_trm.sort.bam
drwxrwxr-x 3 msq msq       4096 Mar  3 11:06 index
-rw-rw-r-- 1 msq msq 1718559500 Mar  4 10:53 P.final.bam
-rw-rw-r-- 1 msq msq 6520916926 Mar  3 17:42 P_trm.bam
-rw-rw-r-- 1 msq msq 4710012734 Mar  3 21:23 P_trm.sort.bam

可参考：对CHIP-seq数据call peaks应该选取unique比对的reads吗？

http://www.bio-info-trainee.com/1869.html

用MACS2软件call peaks

可以call narrow peak(默认）和broad peak比较一下

#随便起的文件名reg，我们做的单样本，不用control
macs2 callpeak -t Example.bed -g mm --nomodel --shift -100 --extsize 200 -n reg --outdir ../peak/narrowpeak/

#broad peak calling
macs2 callpeak -t Example.bed -g mm --nomodel --shift -100 --extsize 200 --broad --broad-cutoff 0.1 -n broad --outdir ../peak/broadpeak/

#批量shell脚本，如果没创建peaks文件夹会自动创建
ls *.bed | while read id ;do (macs2 callpeak -t $id  -g mm --nomodel --shift  -100 --extsize 200  -n ${id%%.*} --outdir ../peaks/) ;done 
#如果不确定代码对不对可以加个echo先不运行：
ls *.bed | while read id ;do (echo macs2 callpeak -t $id  -g mm --nomodel --shift  -100 --extsize 200  -n ${id%%.*} --outdir ../peaks/) ;done 
#

具体参数见：https://www.jianshu.com/p/53d97099b739

和：https://www.jianshu.com/p/21e8c51fca23

同样在开头会显示参数信息

INFO  @ Wed, 04 Mar 2020 11:04:32: 
# Command line: callpeak -t D.bed -g mm --nomodel --shift -100 --extsize 200 -n D --outdir ../peaks/
# ARGUMENTS LIST:
# name = D
# format = AUTO
# ChIP-seq file = ['D.bed']
# control file = None
# effective genome size = 1.87e+09
# band width = 300
# model fold = [5, 50]
# qvalue cutoff = 5.00e-02
# The maximum gap between significant sites is assigned as the read length/tag size.
# The minimum length of peaks is assigned as the predicted fragment length "d".
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 10000 bps
# Broad region calling is off
# Paired-End mode is off

broad peak和narrow peak输出文件会有区别

$ ll
total 15276
-rw-rw-r-- 1 msq msq 1925114 Mar  4 11:25 D_peaks.broadPeak
-rw-rw-r-- 1 msq msq 2062548 Mar  4 11:25 D_peaks_broad.xls
-rw-rw-r-- 1 msq msq 2966429 Mar  4 11:25 D_peaks.gappedPeak
-rw-rw-r-- 1 msq msq 2410497 Mar  4 11:29 P_peaks.broadPeak
-rw-rw-r-- 1 msq msq 2579887 Mar  4 11:29 P_peaks_broad.xls
-rw-rw-r-- 1 msq msq 3686955 Mar  4 11:29 P_peaks.gappedPeak
$ ll
total 12704
-rw-rw-r-- 1 msq msq 2046214 Mar  4 11:08 D_peaks.narrowPeak
-rw-rw-r-- 1 msq msq 2332313 Mar  4 11:08 D_peaks.xls
-rw-rw-r-- 1 msq msq 1289284 Mar  4 11:08 D_summits.bed
-rw-rw-r-- 1 msq msq 2646609 Mar  4 11:11 P_peaks.narrowPeak
-rw-rw-r-- 1 msq msq 3012054 Mar  4 11:11 P_peaks.xls
-rw-rw-r-- 1 msq msq 1667105 Mar  4 11:11 P_summits.bed

deeptools可视化

首先把bam文件转成bw文件，载入igv可视化。igv先打开载入mm10参考基因组

​ bamCoverage 利用测序数据比对结果转换为基因组区域reads覆盖度结果。可以自行设定覆盖度计算的窗口大小(bin)；bamCoverage 内置了各种标准化方法：scaling factor, Reads Per Kilobase per Million mapped reads (RPKM), counts per million (CPM), bins per million mapped reads (BPM) and 1x depth (reads per genome coverage, RPGC).

链接：https://www.jianshu.com/p/2b386dd437d3

cd  ~/project/atac/new
source activate atac
#ls  *.bam  |xargs -i samtools index {} 
ls *final.bam |while read id;do
nohup bamCoverage -p 5 --normalizeUsing CPM -b $id -o ${id%%.*}.final.bw & 
done 
ls *final.bam |while read id;do (bamCoverage -p 5 --normalizeUsing CPM -b $id -o ${id%%.*}.final.bw);done

查看TSS附件信号强度：

注意：computeMatrix给基因组区段打分，产生的文件可用于plotHeatmap 和plotProfiles作图；基因组区段可以是基因或其他区域，使用BED格式文件定义即可（需要给一个参考的注释文件，也就是bed文件，基本的只需要三列：chrom(染色体或contig），chromStart，chromEnd(对应特征的起始位置）三列必选）。

特别注意：bed文件坐标为一半开半闭区间[start, end),所以如果是[10,20)，实际上只提取了10,11，…19 这十个位点，对应ucsc上的即为染色体坐标的10-19位碱基。ucsc上染色体坐标也是从0开始。

自己去ucsc下载对应区段：https://genome.ucsc.edu/cgi-bin/hgTables

点击Get output之后会给我们一个选择输出形式的对话框，在Create one BED record per下面有一些选项，比如这里默认是Whole Gene，当然我们也可以选择启动子区域、外显子加周边区域、5' UTR区域、3' UTR区域等生成我们想要的BED文件。我选的默认

具体参考deeptools使用指南：https://www.jianshu.com/p/2b386dd437d3

## both -R and -S can accept multiple files 
mkdir -p  ~/project/atac/tss
cd   ~/project/atac/tss 
source activate atac
computeMatrix reference-point  --referencePoint TSS  -p 5  \
-b 10000 -a 10000    \   #这里表示画上下游10k
-R ~/project/atac/tss/ucsc.refseq.bed  \
-S ~/project/atac/reference/*.bw  \
--skipZeros  -o matrix1_test_TSS.gz  \
--outFileSortedRegions regions1_test_genes.bed

##     both plotHeatmap and plotProfile will use the output from   computeMatrix
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720 

### 如果要批处理 ，需要学习好linux命令。

前方R报错高能

另外提前准备在浏览器的Rstudio里下载chipseeker的包时，报错了：

g++: internal compiler error: Killed (program cc1plus)
Please submit a full bug report,
with preprocessed source if appropriate.
See <http://bugzilla.redhat.com/bugzilla> for instructions.
make: *** [fastLm.o] Error 4
ERROR: compilation failed for package ‘RcppEigen’

显示依赖包编译失败，考虑到之前我把最开始配置的软件很多为了省存储就删掉了，现在只能重新再安装一遍看能不能解决问题。

然而现在有一个新的问题：

Warning in install.packages(...) :
  'lib = "/usr/lib64/R/library"' is not writable

？？？权限问题？于是我打算以管理员身份重新登录R操作，然而。。。。我发现root根本登录不上Rstudio，本来怀疑密码记错，但是在终端切换是可以的，所以？这个问题我现在还没解决

主要是下面这个，以为换镜像就可以了，但是还是不行，我测试了其他包时可以下载的，只能试从github上下载,从git上下载包要先下载devtools，

Loading required package: ChIPseeker
Bioconductor version 3.10 (BiocManager 1.30.10), R 3.6.0 (2019-04-26)
Installing package(s) 'BiocVersion', 'ChIPseeker'
Error in readRDS(dest) : error reading from connection#没有找到解决方案

install.packages("devtools")
install_github(" YuLab-SMU/ChIPseeker")#报错

尝试所有方案最后还是报错了，看来只能用自己本地的R做

不过最后我找到了问题所在：

#修改了.libPath()后我发现所有的安装错误都是这个包引起的，报错中的ERROR: compilation failed for package ‘RcppEigen’
#于是我重新指定了libpath,且指定安装位置如下，但还是安装失败
install.packages("RcppEigen",lib = "/home/msq/R/x86_64-redhat-linux-gnu-library/3.6")
#ubuntu上是可以用apt-get直接下载R包的，但是CentOs好像不行，我能想到的最终解决方案是：直接把本地的这个包传上去，但是传上去以后,library(RcppEigen)还是报错的
#当我又想放弃时，突然想到其实我的重点不是在这个包，或许我只需要让R认为这个包存在就好了
system call failed: Cannot allocate memory
#上面这个报错不止一次出现，可能是RAM不够造成的，毕竟才2G
#真的成功了！
library(ChIPseeker)

Registered S3 method overwritten by 'enrichplot':
  method               from
  fortify.enrichResult DOSE
ChIPseeker v1.22.1  For help: https://guangchuangyu.github.io/software/ChIPseeker

If you use ChIPseeker in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Qing-Yu He. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 2015, 31(14):2382-2383
#得出结论：
it looks like R doesn't get along very well with low RAM machines

推荐看看这篇linux对memmory和swap的讲解：https://www.cnblogs.com/NBYSY/p/11315902.html

总结教训就是

• 对于报错一定要认真思考和查看帮助，CRAN包官网上会有一些issues供参考。
• R语言对内存要求比想象的高

在linux下载软件的相互依赖非常头疼（conda每次下载一个软件甚至可以给你一口气下几个G，然后你就没空间放数据了（土豪除外））

所以不想耗太多时间在这方面的新手建议直接用本机R做下游

后台运行computeMatrix的情况，运行这个命令时经常服务器中断，可能还是服务器性能不足的问题

3935 msq       20   0  380344 214544  10496 R 98.3 11.4   4:48.33 computeMatrix       

#看genebody区的信号强度
computeMatrix scale-regions  -p 15 -R ~/project/atac/tss/ucsc.refseq.bed -S ~/project/atac/reference/*.final.bw -b 2000 -a 2000 --skipZeros -o matrix1_test_body.gz             

peaks注释

因为上步的compute小服务器的计算力好像做不了，不影响下游所以先跳过，拿到之前做的narrow和broad peaks，放在R里面先进行注释

另外在颉伟老师的paper有很多peaks文件供下载: GSE124718

R运行anno时的截图：feature的概念和gene相当，从Seurat2到Seurat3的变迁可以看出

preparing features information...        2020-03-05 14:03:56 
identifying nearest features...      2020-03-05 14:03:59 
calculating distance from peak to TSS...     2020-03-05 14:04:00 
assigning genomic annotation...      2020-03-05 14:04:00 
adding gene annotation...            2020-03-05 14:04:30 
Loading required package: org.Mm.eg.db

'select()' returned 1:many mapping between keys and columns

assigning chromosome lengths             2020-03-05 14:04:31 
done...                     2020-03-05 14:04:31 

需要很好地理解bioconductor很多包内置函数创建的对象，对象可以理解成一个高维的数据框，对数据框的操作同样可以很好地应用到对象

统计peak在promoter，exon，intron和intergenic区域的分布

rm(list = ls())

bedPeaksFile         = '../narrow/D_summits.bed'; 
bedPeaksFile
## loading packages
require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile ) 
#排除chr_random和chrUn
keepChr= !grepl('_',seqlevels(peak))
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Mm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter) 
# 然后查看这些peaks在所有基因的启动子附近的分布情况，热图模式
heatmap<- tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
ggsave(heatmap,file = 'heatmap.pdf')
# 然后查看这些peaks在所有基因的启动子附近的分布情况，信号强度曲线图，饼图
signal<- plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
ggsave(signal,file = "signal.pdf")
plotAnnoPie(peakAnno)

heatmap效果比较差，但我重新做broad peak注释时发现结果比较明显，每一行有相对宽的红色peaks

另外，broad peaks的信号强度曲线图明显平滑，相比起narrow peak尖锐的峰

table(tagMatrix)#可以看到大部分都是0值，稀疏矩阵
tagMatrix
       0        1 
67419342    13895

上述anno代码运行结果写成了Rmd报表，可在报表里查看。

用Y叔的ChIPseeker对peaks进行注释和可视化：https://www.jianshu.com/p/26aaba19a605

Motif 寻找及注释

Homer

Homer软件的介绍-最全面而详细的找motif教程：https://mp.weixin.qq.com/s/U7iAS75Mlg6aJFqGp1LfDg

Homer依赖perl，这个软件最初安装非常麻烦，有conda会好很多，但仍然需要进行配置：注意要找的配置文件在share下，不要用which homer找

perl ~/miniconda3/envs/atac/share/homer-4.10-0/configureHomer.pl -install mm10 
#下载可能会比较慢这个看运气了，下载好后，查看一下已有的数据包：
perl ~/miniconda3/envs/atac/share/homer-4.10-0/configureHomer.pl -list
#截了部分结果，含有organism的ontology,promoter，genomes and annotation information，而且物种很丰富
PROMOTERS
-   yeast-p v5.5    yeast promoters (yeast)
-   mouse-p v5.5    mouse promoters (mouse)
-   zebrafish-p v5.5    zebrafish promoters (zebrafish)
-   rat-p   v5.5    rat promoters (rat)
-   arabidopsis-p   v6.3    arabidopsis promoters (arabidopsis)
-   frog-p  v5.5    frog promoters (frog)
-   fly-p   v5.5    fly promoters (fly)
-   human-p v5.5    human promoters (human)
-   worm-p  v5.5    worm promoters (worm)
-   chicken-p   v5.5    chicken promoters (chicken)

#install的参考基因组放在genomes文件夹里
ls -lh  ~/miniconda3/envs/atac/share/homer-4.9.1-5/data/genomes  
mkdir -p  ~/project/atac/motif
cd   ~/project/atac/motif  
source activate atac#有的conda激活环境的命令是conda activate 'env'
ls ../peaks/*.narrowPeak |while read id;
do 
file=$(basename $id )
sample=${file%%.*} 
echo $sample #awk命令随心所欲提取对应的列，找motif理论上只需要对应peaks的起始位置和编号即可
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id > ${sample}.homer_peaks.tmp
nohup findMotifsGenome.pl ${sample}.homer_peaks.tmp  mm10 ${sample}_motifDir -len 8,10,12  & 
done 
#找motif依然需要相对大的计算资源，而且运行后开头的info会自动检查你的peaks质量，因为就两个样本没有用nohup
Position file = D_peaks.homer_peaks.tmp
    Genome = mm10
    Output Directory = homer_motifDir
    Motif length set at 8,10,12,
    Found mset for "mouse", will check against vertebrates motifs
    Peak/BED file conversion summary:
        BED/Header formatted lines: 0
        peakfile formatted lines: 28430

    Peak File Statistics:
        Total Peaks: 28430
        Redundant Peak IDs: 0
        Peaks lacking information: 0 (need at least 5 columns per peak)
        Peaks with misformatted coordinates: 0 (should be integer)
        Peaks with misformatted strand: 0 (should be either +/- or 0/1)

    Peak file looks good!
....#省略一堆信息
Reading input files...
    56860 total sequences read
    Cache length = 11180
    Using binomial scoring  #打分模型
#感觉跑的时候很费劲毕竟内存只有那么大

以下是我的输出文件：可以打开html报表查看详细信息，我把输出文件都转移到了本地电脑对应proj文件夹下方便查看

(atac) msq 11:41:58 ~/project/atac/motif/homer_motifDir
$ ll
total 1152
-rw-rw-r-- 1 msq msq  18826 Mar  8 11:16 homerMotifs.all.motifs
-rw-rw-r-- 1 msq msq  10043 Mar  8 10:46 homerMotifs.motifs10
-rw-rw-r-- 1 msq msq      0 Mar  8 10:46 homerMotifs.motifs12
-rw-rw-r-- 1 msq msq   8783 Mar  8 10:32 homerMotifs.motifs8
drwxrwxr-x 2 msq msq   4096 Mar  8 11:19 homerResults
-rw-rw-r-- 1 msq msq  79851 Mar  8 11:19 homerResults.html
drwxrwxr-x 2 msq msq  20480 Mar  8 10:30 knownResults
-rw-rw-r-- 1 msq msq 970164 Mar  8 10:30 knownResults.html
-rw-rw-r-- 1 msq msq  48015 Mar  8 10:30 knownResults.txt
-rw-rw-r-- 1 msq msq     63 Mar  8 10:11 motifFindingParameters.txt
-rw-rw-r-- 1 msq msq   1888 Mar  8 10:14 seq.autonorm.tsv

在线网站工具

如果觉得太耗时，打算用其他方法做，参考：

meme 在线网站（有丰富的套件工具，enrich推荐CentriMo，比如我对broad peak的出图效果，还可导出tsv文件:http://meme-suite.org/opal-jobs/appCENTRIMO_5.1.11583556586142956512218/centrimo.html#data_sec）

meme suite官网 ：http://meme-suite.org/tools/meme，需要 获取 序列参考代码，参考下面这个R脚本(需要下载对应的BSgenome的包，还是看网速)

https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step7-peaks2sequence.R

R package downstream analysis

注意getAllpeakssequence函数需要对应的BSgenome对象，如果用ChIPpeakAnno做注释，可能要借助其他GRanges对象来作为AnnotationData, 在包的下载位置，如TxDb.Mmusculus.UCSC.mm10.knownGene package，可以在网页文档介绍中找到包的数据来源和update时间）

丰富高能的R包，比如 motifmatchr包 也可以做，加粗表示下游分析还有更多, 比如ChIPpeakAnno还可以找enhancer:Find possible enhancers with DNA interaction data，DNA 的interactive block可以达到100k；还有comparing binding profiles from multiple transcription factors (TFs)：Given two or more peak lists from different TFs, one may be interested in finding whether DNA binding profile of those TFs are correlated, and if correlated, what is the common binding pattern. The workflow here shows how to test the correlation of binding profiles of three TFs and how to discover the common binding pattern. also Heatmaps showing Tf binding comparison

motifmatchr： https://bioconductor.org/packages/release/bioc/html/motifmatchr.html

MotIV:http://bioconductor.org/packages/devel/bioc/html/MotIV.html

rGADEM:http://bioconductor.org/packages/devel/bioc/html/rGADEM.html

ChIPpeakAnno(集成下游分析平台，含示例）：http://www.bioconductor.org/packages/release/bioc/vignettes/ChIPpeakAnno/inst/doc/pipeline.html

学习一下ChIPpeakAnno(一)：https://www.baidu.com/link?url=Tn_wFzlQIARFh4HePlMdpK2i29m6KhzBPHDr-GKdQkshwk3IqIy5snrS1qbcI1mf&wd=&eqid=e669443900010707000000065e635dcc

GREAT: Genomic Regions Enrichment of Annotations Tool：http://great.stanford.edu/public/html/splash.php regulatory domains identification and GO annotation

rGREAT:http://www.bioconductor.org/packages/release/bioc/vignettes/rGREAT/inst/doc/rGREAT.html

参考：自学CHIP-seq分析第七讲~peaks注释：http://www.bio-info-trainee.com/1752.html

ref

可以参考这篇文章的方法：

paper:ChIP-Enrich: gene set enrichment testing for ChIP-seq data:https://academic.oup.com/nar/article/42/13/e105/1275211

拿到bam文件，calling peaks文件和分组信息后可以尝试这个包：

第5篇：对ATAC-Seq/ChIP-seq的质量评估（二）——ChIPQC:https://www.jianshu.com/p/a11147808d14

生信技能树以前的推文：一篇文章学会ChIP-seq分析：https://mp.weixin.qq.com/s/tnbhwZ16QiAnWulCvBNFsQ

关于如何鉴定增强子：

• https://www.cnblogs.com/leezx/p/10812942.html
• http://www.360doc.com/content/18/1216/11/52645714_802150818.shtml

如何对基因组序列进行注释：https://www.jianshu.com/p/931e9821c45a