寫在前面

可能第一眼看到這個標題大家都會很奇怪蔑匣，因為2G內(nèi)存的服務(wù)器還要分析上游的話咽块，不太可能潮剪，確實在中間我把比對的過程跳過了
服務(wù)器生信分析平臺配置參考https://zhuanlan.zhihu.com/p/74423987
生信是實踐性很強的交叉學(xué)科，通過實踐提高自己思考和解決問題的能力

本文結(jié)合linux操作和R叹洲，R單獨在本地其實也可以做柠硕，用esATAC，一步包裝到位就是

主要參考文章是生信技能樹上的教程：http://www.reibang.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里面比較簡單，上下游結(jié)合的教程整理蝗柔；
• TCGA 數(shù)據(jù)庫挖掘闻葵，推薦在cellpress上的泛癌研究文獻：https://www.cell.com/pb-assets/consortium/pancanceratlas/pancani3/index.html
• GEO-master的airwayRNA
• ChIPpeakAnno, atacseq proj分析報告
• 加拿大生物信息學(xué)研討會資源寶藏：https://bioinformatics.ca/workshops/2018-epigenomic-data-analysis/
  • 對應(yīng)youtube鏈接：https://www.youtube.com/channel/UCKbkfKk65PZyRCzUwXOJung
• 單細胞培訓(xùn)：https://www.csc.fi/web/training/-/scrnaseq可以下載code

在win10上用linux子系統(tǒng)做RNAseq

參考的視頻：https://www.bilibili.com/video/av65970178?p=4

參考筆記和paper: http://www.reibang.com/p/6fa6da3bae40

膽子夠大可以嘗試搭建linux虛擬機

下面是幾個points：

1. win10默認掛載到/mnt上，這對于拷貝文件很方便癣丧，不用FileZilla傳來傳去

2. 下載的軟件zip槽畔、安裝的軟件路徑、編譯軟件等目錄要安排好,養(yǎng)成好習(xí)慣胁编，conda配置環(huán)境work時也是如此:

   *bin：可執(zhí)行程序厢钧，添加所有軟件鏈接方便加入環(huán)境變量*
   
   *biosoft：軟件安裝目錄，解壓編譯后放這里*
   
   *src：軟件源代碼嬉橙，最開始下載的包放在這里*
   

   下載軟件失敗早直，發(fā)現(xiàn)原來在備份sources.list的時候把sources打成了source

3. 操作.bashrc文件重新個性化配置時（用戶的.bashrc里放著用戶自己的環(huán)境變量,export $PATH可以添加軟件路徑到里面，推薦用conda解決大部分軟件環(huán)境問題）市框，vim的操作要熟悉

4. 編譯gffcompare源代碼的時候霞扬，g++編譯器下載失敗，apt-get做不了枫振，說是有一系列依賴包沒下喻圃，然后我發(fā)現(xiàn)那些依賴包也下不了，想放棄的時候從頭檢查感覺可能是etc/apt下的sources.list文件有問題粪滤，把阿里云的軟件下載鏡像站點都換成清華的以后斧拍，編譯gffcompare就成功了

5. HISAT2速度相當(dāng)快,stringtie重構(gòu)轉(zhuǎn)錄本同時組裝轉(zhuǎn)錄本和評估其表達水平，gffcompare評估組裝效果（檢查拼接的轉(zhuǎn)錄本有多少匹配到了參考基因組注釋文件杖小，gffcompare可以用來compare, merge, annotate and estimate accuracy of one or more GFF files), 最后再Stringtie估算轉(zhuǎn)錄本表達量并生成Ballgown能使用的統(tǒng)計表格肆汹，后面就是下游的事了

   
   
   stringtie.png

如上圖，stringtie是新轉(zhuǎn)錄本窍侧，可能是lncRNA

偶然間看到騰訊云特價的學(xué)生用服務(wù)器后果斷買了县踢，最便宜的1核2G學(xué)生用服務(wù)器，盡管遠遠不夠生信分析的要求伟件，體驗還是可以的，有條件的同學(xué)建議買16核64G服務(wù)器议经，如果想跑很多個樣本的數(shù)據(jù)的話（比如單細胞）斧账。

下載數(shù)據(jù)折騰3-4天

因為網(wǎng)速/有的站點網(wǎng)絡(luò)不穩(wěn)定的問題，下載數(shù)據(jù)很多時候也很麻煩

國內(nèi)站點下載數(shù)據(jù)一直是個問題煞肾，后來才看到有人推薦用海外云服務(wù)器咧织，按用量計費（可能比較劃算），手動下載是不現(xiàn)實的（不支持用wget籍救，用axel多線程斷點續(xù)傳花了兩天才下了一個3G的測試數(shù)據(jù)（HISAT文件的RNAseq pipeline里的）)聽說aspera超高速下載sra數(shù)據(jù)很棒习绢，嘗試了一下失敗（ncbi可能取消了ascp協(xié)議，或者我自己配置的問題闪萄，不支持ncbi下載梧却，還有一個重要原因是ncbi的ftp目錄組織不太好，經(jīng)常會變化）但是ncbi不是唯一的出路

最后選擇了EBI-ERA數(shù)據(jù)庫下載败去，同樣支持SRR號下載放航，而且界面科學(xué)，點擊ftp復(fù)制鏈接即可圆裕，而且下載后直接是fasta文件广鳍，省去了ncbi下載的sra文件還要fastq-dump的步驟，而且ERA是絕對支持ascp的吓妆！速度超快幾百M/s赊时，即使數(shù)據(jù)量很大但也沒花多少時間，這個時候考慮到我的硬盤要省著點用（因為后面測試時因存儲問題行拢，文件是truncated的祖秒，只好刪掉兩個測試數(shù)據(jù)）

#注意密鑰文件位置，vol1后根據(jù)自己要下載數(shù)據(jù)的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

下載數(shù)據(jù)在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581

conda配置環(huán)境

照搬代碼剂陡，不解釋狈涮，注意是py2，安裝完一定要檢查軟件有沒有安裝成功鸭栖，否則像我一樣后面又刪除環(huán)境重裝歌馍。添加channel后看一下channel的優(yōu)先級，建議把bioconda放最上面晕鹊，減少search時間松却。

關(guān)于conda的基礎(chǔ)操作：http://www.reibang.com/p/edaa744ea47d

推薦安裝軟件的課程：https://www.bilibili.com/video/av82983858?p=7

所有的軟件安裝失敗都是環(huán)境配置造成的

# 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后需要設(shè)置鏡像。
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這個環(huán)境下面
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文件，否則時間長了特別占存儲

質(zhì)控

fastqc + multiqc兩件套飞几，不多解釋砚哆，可以好好看看輸出的報告文件，沒安裝的自己安裝

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

fastqc官網(wǎng)屑墨，含示例報告文件：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

推薦trim_galore和shell編程

trim_galore去接頭躁锁，這個軟件非常智能，自動探測adaptor，比trimmomatic好

參考該軟件：http://www.reibang.com/p/7a3de6b8e503

參考優(yōu)秀的shell筆記：http://www.reibang.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

當(dāng)時運行的界面開頭是一些運行參數(shù)，可以看到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最好在當(dāng)前都是fastqc report的路徑下
multiqc .

比對

比對計算量相當(dāng)大，這里我的服務(wù)器超內(nèi)存了：processing memory不足以運行bowtie2（忘記截屏）

wget下載索引還是比較慢，最好不要自己構(gòu)建index

# 索引大小為3.2GB刁标， 不建議自己下載基因組構(gòu)建颠通，可以直接下載索引文件，代碼如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip

因為我自己服務(wù)器沒法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到后#面我發(fā)現(xiàn)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分數(shù)

[SAMtools] 常用指令總結(jié)

理解并操作bam文件：

很多NGS文件都分成header和record部分，header里面還會包括此文件生成的過程信息,沒有sort的bam文件在開頭是這樣的：

@HD VN:1.0 SO:unsorted

sort后：

@HD VN:1.0 SO:coordinate

header內(nèi)容不多匈织，也不會太復(fù)雜浪默，每一行都用‘@’ 符號開頭，里面主要包含了版本信息缀匕，序列比對的參考序列信息纳决，如果是標準工具（bwa，bowtie乡小，picard）生成的BAM阔加，一般還會包含生成該份文件的參數(shù)信息（如上圖），@PG標簽開頭的那些满钟。這里需要重點提一下的是header中的@RG也就是Read group信息胜榔，這是在做后續(xù)數(shù)據(jù)分析時專門用于區(qū)分不同樣本的重要信息。它的重要性還體現(xiàn)在湃番，如果原來樣本的測序深度比較深夭织，一般會按照不同的lane分開比對，最后再合并在一起吠撮，那么這個時候你會在這個BAM文件中看到有多個RG尊惰，里面記錄了不同的lane，甚至測序文庫的信息泥兰，唯一不變的一定是SM的sample信息择浊，這樣合并后才能正確處理。

所以如果一個處理組有好幾個重復(fù)逾条，還需要把這些bam文件merge一下：samtools merge一個output三個input

寫循環(huán)除了while read id do還可以用xargs命令：

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

參考鏈接：http://www.reibang.com/p/364e640d3c9f

可以對sort后的文件構(gòu)建index，方便查看的時候直接跳到自己感興趣的區(qū)域查看：

samtools index in.bam  # 生成in.bam的索引文件in.bam.bai
samtools view in.bam chr22            # 跳轉(zhuǎn)到chr22染色體
samtools view in.bam chr22:16050103   # 跳轉(zhuǎn)到chr22:16050103位置
samtools view in.bam chr22:16050103-16050103  # 只查看該位置

另外發(fā)現(xiàn)有一些比較奇特的染色體名：（我用的是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

查教程發(fā)現(xiàn)：

“chr*_random sequences” 知道來自哪條染色體但不知道具體位置的序列

“chrUn_* sequences” 知道來自人類基因組序列投剥，但不知道與染色體的關(guān)系

詳細參考：http://www.reibang.com/p/2cd5e6e310c9

還可以對bam/sam文件做統(tǒng)計：

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

輸出結(jié)果：

##總的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之間的距離符合設(shè)置的閾值
35840362 + 0 properly paired (57.46% : N/A)
#PE中兩條都比對到參考序列上的reads數(shù)
47374920 + 0 with itself and mate mapped
#單獨一條匹配到參考序列上的reads數(shù)，和上一個相加，則是總的匹配上的reads數(shù)
1765705 + 0 singletons (2.83% : N/A)
#PE中兩條分別比對到兩條不同的參考序列的reads數(shù)
1226554 + 0 with mate mapped to a different chr
#PE中兩條分別比對到兩條不同的參考序列的reads數(shù)（比對質(zhì)量>=5）
19944 + 0 with mate mapped to a different chr (mapQ>=5)

過濾BAM文件

ATAC-Seq與其他方法不同的一點是需要過濾去除線粒體（如果是植物吃警，還需要過濾葉綠體）糕篇，因為線粒體DNA是裸露的，也可以被Tn5酶識別切割

參考：http://www.reibang.com/p/5bddf9c9ca4c

過濾唯一比對的reads

#這句話-h指定打印頭文件信息酌心，-t指定線程拌消，-f指定bam文件格式，-F指定過濾條件安券，[XS]是sam的可選域墩崩，描述序列多重比對上時第二個最有比對的得分，只有當(dāng)序列比對到多個位置時才會出現(xiàn)侯勉，因此[XS] == null除去了多重比對的序列鹦筹，not unmapped 和not duplicate分別用來除去未比對上序列和重復(fù)序列
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重復(fù)

這個我沒有做址貌，因為上面數(shù)據(jù)里沒有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參數(shù)反選
#此時就得到了唯一比對且已經(jīng)去除過線粒體的比對文件铐拐，可以用于接下來的peaks calling

將bam文件轉(zhuǎn)為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數(shù)據(jù)call peaks應(yīng)該選取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腳本螟凭，如果沒創(chuàng)建peaks文件夾會自動創(chuàng)建
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 
#

具體參數(shù)見：http://www.reibang.com/p/53d97099b739

和：http://www.reibang.com/p/21e8c51fca23

同樣在開頭會顯示參數(shù)信息

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輸出文件會有區(qū)別

$ 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文件轉(zhuǎn)成bw文件虚青，載入igv可視化。igv先打開載入mm10參考基因組

? bamCoverage 利用測序數(shù)據(jù)比對結(jié)果轉(zhuǎn)換為基因組區(qū)域reads覆蓋度結(jié)果赂摆⌒荆可以自行設(shè)定覆蓋度計算的窗口大小(bin)；bamCoverage 內(nèi)置了各種標準化方法：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).

鏈接：http://www.reibang.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給基因組區(qū)段打分烟号，產(chǎn)生的文件可用于plotHeatmap 和plotProfiles作圖绊谭；基因組區(qū)段可以是基因或其他區(qū)域，使用BED格式文件定義即可（需要給一個參考的注釋文件汪拥，也就是bed文件达传，基本的只需要三列：chrom(染色體或contig），chromStart迫筑，chromEnd(對應(yīng)特征的起始位置）三列必選）宪赶。

特別注意：bed文件坐標為一半開半閉區(qū)間[start, end),所以如果是[10,20)，實際上只提取了10,11脯燃，…19 這十個位點搂妻，對應(yīng)ucsc上的即為染色體坐標的10-19位堿基。ucsc上染色體坐標也是從0開始辕棚。

自己去ucsc下載對應(yīng)區(qū)段：https://genome.ucsc.edu/cgi-bin/hgTables

點擊Get output之后會給我們一個選擇輸出形式的對話框欲主，在Create one BED record per下面有一些選項邓厕，比如這里默認是Whole Gene，當(dāng)然我們也可以選擇啟動子區(qū)域扁瓢、外顯子加周邊區(qū)域详恼、5' UTR區(qū)域、3' UTR區(qū)域等生成我們想要的BED文件引几。我選的默認

具體參考deeptools使用指南：http://www.reibang.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 

### 如果要批處理 昧互，需要學(xué)習(xí)好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’

顯示依賴包編譯失敗敞掘，考慮到之前我把最開始配置的軟件很多為了省存儲就刪掉了，現(xiàn)在只能重新再安裝一遍看能不能解決問題贿讹。

然而現(xiàn)在有一個新的問題：

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

渐逃？？民褂？權(quán)限問題茄菊？于是我打算以管理員身份重新登錄R操作，然而赊堪。面殖。。哭廉。我發(fā)現(xiàn)root根本登錄不上Rstudio脊僚，本來懷疑密碼記錯，但是在終端切換是可以的遵绰，所以辽幌？這個問題我現(xiàn)在還沒解決

主要是下面這個，以為換鏡像就可以了椿访，但是還是不行乌企，我測試了其他包時可以下載的，只能試從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()后我發(fā)現(xiàn)所有的安裝錯誤都是這個包引起的，報錯中的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)還是報錯的
#當(dāng)我又想放棄時陋葡，突然想到其實我的重點不是在這個包，或許我只需要讓R認為這個包存在就好了
system call failed: Cannot allocate memory
#上面這個報錯不止一次出現(xiàn)彻采，可能是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
#得出結(jié)論：
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

總結(jié)教訓(xùn)就是

• 對于報錯一定要認真思考和查看幫助，CRAN包官網(wǎng)上會有一些issues供參考柴梆。
• R語言對內(nèi)存要求比想象的高

在linux下載軟件的相互依賴非常頭疼（conda每次下載一個軟件甚至可以給你一口氣下幾個G，然后你就沒空間放數(shù)據(jù)了（土豪除外））

所以不想耗太多時間在這方面的新手建議直接用本機R做下游

后臺運行computeMatrix的情況终惑，運行這個命令時經(jīng)常服務(wù)器中斷绍在，可能還是服務(wù)器性能不足的問題

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

#看genebody區(qū)的信號強度
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小服務(wù)器的計算力好像做不了，不影響下游所以先跳過雹有，拿到之前做的narrow和broad peaks偿渡，放在R里面先進行注釋

另外在頡偉老師的paper有很多peaks文件供下載: GSE124718

R運行anno時的截圖：feature的概念和gene相當(dāng)，從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很多包內(nèi)置函數(shù)創(chuàng)建的對象霸奕，對象可以理解成一個高維的數(shù)據(jù)框溜宽，對數(shù)據(jù)框的操作同樣可以很好地應(yīng)用到對象

統(tǒng)計peak在promoter，exon质帅，intron和intergenic區(qū)域的分布

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注釋時發(fā)現(xiàn)結(jié)果比較明顯，每一行有相對寬的紅色peaks

另外魄揉，broad peaks的信號強度曲線圖明顯平滑剪侮，相比起narrow peak尖銳的峰

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

上述anno代碼運行結(jié)果寫成了Rmd報表洛退，可在報表里查看瓣俯。

用Y叔的ChIPseeker對peaks進行注釋和可視化：http://www.reibang.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 
#下載可能會比較慢這個看運氣了摇零，下載好后推掸，查看一下已有的數(shù)據(jù)包：
perl ~/miniconda3/envs/atac/share/homer-4.10-0/configureHomer.pl -list
#截了部分結(jié)果，含有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激活環(huán)境的命令是conda activate 'env'
ls ../peaks/*.narrowPeak |while read id;
do 
file=$(basename $id )
sample=${file%%.*} 
echo $sample #awk命令隨心所欲提取對應(yīng)的列，找motif理論上只需要對應(yīng)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質(zhì)量毡泻，因為就兩個樣本沒有用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  #打分模型
#感覺跑的時候很費勁畢竟內(nèi)存只有那么大

以下是我的輸出文件：可以打開html報表查看詳細信息，我把輸出文件都轉(zhuǎn)移到了本地電腦對應(yīng)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

在線網(wǎng)站工具

如果覺得太耗時粘优，打算用其他方法做仇味，參考：

meme 在線網(wǎng)站（有豐富的套件工具呻顽，enrich推薦CentriMo，比如我對broad peak的出圖效果丹墨，還可導(dǎo)出tsv文件:http://meme-suite.org/opal-jobs/appCENTRIMO_5.1.11583556586142956512218/centrimo.html#data_sec）

meme suite官網(wǎng) ：http://meme-suite.org/tools/meme廊遍，需要 獲取 序列參考代碼，參考下面這個R腳本(需要下載對應(yīng)的BSgenome的包贩挣，還是看網(wǎng)速)

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

R package downstream analysis

注意getAllpeakssequence函數(shù)需要對應(yīng)的BSgenome對象喉前，如果用ChIPpeakAnno做注釋，可能要借助其他GRanges對象來作為AnnotationData, 在包的下載位置王财，如TxDb.Mmusculus.UCSC.mm10.knownGene package卵迂，可以在網(wǎng)頁文檔介紹中找到包的數(shù)據(jù)來源和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

學(xué)習(xí)一下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

參考：自學(xué)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的質(zhì)量評估（二）——ChIPQC:http://www.reibang.com/p/a11147808d14

生信技能樹以前的推文：一篇文章學(xué)會ChIP-seq分析：https://mp.weixin.qq.com/s/tnbhwZ16QiAnWulCvBNFsQ

關(guān)于如何鑒定增強子：

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

如何對基因組序列進行注釋：http://www.reibang.com/p/931e9821c45a