Ⅰ chapter1.高通量序列實驗簡介：設計與生物信息學分析
2021/4/16
1、設計高通量測序實驗的步驟
2、介紹了最廣泛使用的應用，並描述了基本的測序概念。
3、可用於生物信息學分析的各種軟體程序，以理解測序數據。
1、Insert ：用於測序的DNA片段
2、Read ： Insert 被測序到的部分
3、Single Read（SR） ：一種只從 Insert 序列一端測序的測序程序
4、Pair Read（PR） ：一種從 Insert 序列兩端測序的測序程序
5、Flowcell ：連接DNA晶元並進行測序的一種小玻璃晶元。 Flowcell 被探針覆蓋，允許與DNA片段連接的接頭雜交。
6、Lane ： Flowcell 由8個物理分離的通道組成，稱為 Lane 。在所有 Lane 上並行進行測序。
7、Multiplexing/Demultiplexing ：在同一 Lane 上對幾個茄燃樣本進行測序稱為多路復用 Multiplexing ，在一條 Lane 上測序的 Reads 的分離稱為分路復用 Demultiplexing ，通過一個識別每個Reads*的索引將其與已知樣本的索引進行比較。
8、Pipeline ：一系列的計算過程
（一）reading
1、Resequencing ：在一個給定的樣本中找到相對於參考基因組的變體
實驗細節 ：從相關細胞中提取DNA，進行由DNA碎片化和測序組成的樣品制備
基本分析總結 ：將序列片段映射到參考基因組，並通過總結片段與其基因組位點的差
異來識別相對於參考基因組的變異對應的「地圖」
2、Target-enriched sequencing ：靶點富集測序是一種特定的 Resequencing 形式，只
關注特定的基因組基因座。
實驗細節 ：在從細胞中提取DNA並進行樣品制備後，進行一個富集過程來捕獲相關的
位點，靶富集可以使用「定製的」靶富集探針在基因組的特定區域進行，或
使用可用的試劑盒，如exome-enrichment kits。
基本分析總結 ：與 Resequencing 相同
3、De novo assembly ：識別一個基因組序列，而無需任何額外的參考
實驗細節 ：與 Resequencing 相同
基本分析總結 ：組裝過程依賴於DNA片段的重疊。這些重疊被合並成一致序列，稱為
contigs 和 scaffolds 。
（二）counting
1、ChIP-Seq/RIP-Seq ：找到RNA或DNA結合蛋白的結合位置
實驗細節 ：（1）首先，進行了ChIP/RIP實驗：蛋白質與DNA/RNA結合，並與之交
聯。然後DNA/RNA被分裂。
（2）蛋白質 pull down 經歷免疫沉澱過程，交聯被逆轉
（3）對富集於蛋白結陸納虛合位點中的DNA/RNA片段進行測序
基本分析總結 ：被序列排列的片段被映射到基因組中。基因組中豐富的位置是通過檢
測基因組的映射片段的「 peaks (峰)」發現的,這些峰值應該明顯高於在
周圍的位點中已映射的片段，並且與對照樣本相比要高得多----通常
是ChIP實驗的輸入DNA或其他由非特異性抗體進行的免疫沉澱樣
本。
2、RNA-Seq ：檢測和比較基因表達水平
實驗細節 ：從細胞中提取總RNA，在樣品制備過程中，mRNA被 pull down 並破碎。
然後，mRNA片段被逆轉錄成cDNA，cDNA片段測序。
基本分析總結 ：cDNA片段被映射到參考基因組中。映射到每個基因的片段被計數和
標准化，以便比較不同的基因和不同的樣本。在一個RNA-Seq實驗
中，通過檢測映射到一個未注釋區域的基因組上的片段束，可以找到
未標記的基因和轉錄本。
（三）reading/counting
microRNA-Seq ：檢測和計數 microRNAs
實驗細節 ：從細胞中提取總RNA，通過識別大多數已知的microRNA分子共同的自然
早燃 結構來分離microRNA，然後對microRNA片段進行逆轉錄和測序。
基本分析總結 ：被測序的片段被映射到基因組中，然後，微RNA可以被檢測和計數。
1、在 reading 中，覆蓋范圍對應於 平均覆蓋基因組中每個鹼基 的 reads 數量。
一般來說， 30X覆蓋率 被認為是識別基因組變異的最小值，而 de novo 通常需要一個更高的覆蓋范圍。
2、在 counting 中，覆蓋的概念並不簡單，因為the number of reads along the genome is not expected to be uniform.
可幫助評估是否有足夠的reads序列的分析是「 *saturation report* (飽和度報告)」，使用所有的reads確定表達水平，表達水平與 取一部分reads重新計算的表達水平 比較。
1、基因組的重復性
要唯一地對重復區域的read映射進行評分，它必須 比 重復區域 或 邊界相鄰的非重復序列 更長。更長的reads或PE reads允許「拯救」非唯一端，也映射到基因組中的非唯一區域。
2、差異剪接變異
同一基因表達的轉錄本不同時：
3、測序樣本與參考基因組的遺傳距離
如果被測序的樣本與參考基因組有遺傳距離，可能需要更長的reads來確定基因組中每個read的來源。
4、尋找結構變異
基因組的結構變化，如長的插入或缺失，倒位和易位可以通過Paired-End信息找到。
5、De Novo 裝配
挑戰：測序錯誤、低復雜度區域和重復區域等
更長的PE reads會導致更好的裝配，使用一些具有不同insert length的序列庫可以改進組裝過程。
1、 Resequencing：有遺傳距離。。。
2、RNA-Seq：使用來自不同重復的數據，並將其合並為一個具有更高統計顯著性的值。
3、ChIP-Seq：+控制樣本
1、Raw Data 處理
此步驟的可用軟體：Illumina』s CASAVA software，Illumina運行會生成「base-calling」文件(*.bcl)，它們只有在轉換為通用fastq格式時才會非常有用,在此文件轉換過程中，還執行解復用過程，即從同一lane上排序的不同樣本分離讀取。
2、質量控制和read操作
此步驟的可用軟體：CASAVA和FastQC
測序運行完成後，在開始分析之前，應檢查運行的質量是否以下參數，這些參數可能說明樣品和運行的質量。
3、為De Novo Assembly組裝 Contigs 和 Scaffolds
此步驟的可用軟體：SOAPdenovo，ABySS，Velvet，ALL-PATHS
4、mapping
此步驟的可用軟體：BWA ，Bowtie，TopHat
5、 Variant Calling and Filtering
此步驟的可用軟體： SAMtools，GATK，MAQ
幫助檢測變異的兩個基本參數如下：
（1）Coverage at the loci
（2）被測序的等位基因的頻率
6、Assembling Transcripts
7、 Gene Expression Analysis
此步驟的可用軟體：Cufflinks，Myrna
一種常見的歸一化方法FPKM，計算如下：
8、 Peak Detection
此步驟的可用軟體：MACS，SICER
Ⅱ VarBen：NGS變異數據模擬軟體
2021年3月3號，來自 國家衛生健康委臨床檢驗中心、中國科學院計算技術研究所等單位的研究人員開發了腫瘤突變數據模擬軟體VarBen，包括SNV、InDel、CNV、SV等多種變異形式的模擬。同時該軟體支持目前臨床上常用的測序平台，包括Illumina、華大智造MGISEQ以及Ion torrent測序平台。
據介紹，VarBen採用的是比對到參考基因組特定位點的測序reads進行編輯的方式來進行突變模擬，該方法可保留測序過程「濕實驗」部分核酸提取、靶向捕獲、文庫制備以及測序過程中產生的錯誤分布模式，從而保證模擬數據更加的接近真實。研究人員根據不同高通量測序平台的原理以及不同種類基因變異的特點，使得VarBen能夠對幾乎所有的變異類型進行模擬，包括SNV、Indel、Complex insertion-deletion、CNV以及SV。此外，VarBen同時支持全基因組、全外顯子組以及靶向panel測序數據的模擬，並且適用於多個測序平台。
文獻名稱： VarBen: Generating in Silico Reference Data Sets for Clinical Next-Generation Sequencing Bioinformatics Pipeline Evaluation
期刊名稱：The Journal of Molecular Diagnostics （國際分子診斷領域權威期刊） ( IF 5.553 )(JMD)
發表時間：2020-12-18
文章介紹1： https://www.x-mol.com/paper/1340047108494508032
文章介紹2： https://www.sohu.com/a/453971998_120054289
github： https://github.com/nccl-jmli/VarBen
數據路徑： SRP171108
軟體亮點（和其他軟體相比）：
測試軟體：
注意：bwa、samtools、bedtools 要加入環境變數中，流程裡面直接使用的時bwa，samtools、bedtools，沒有給傳軟體路徑的參數
最終生成的bam文件存儲在 outdir/edit.sorted.bam
i)從BAM文件中根據用戶指定的變異等位基因分數(VAF)的指定位置或區域隨機選擇Read；
ii)根據指定的變體類穗旁型(包括)修改、插入或刪除重疊位置或區域的Reads中的SNV和InDel；
iii)編輯後的Reads被重新映射到參考基因組以獲得適當的比對信息並與從原始輸入BAM中剩下的的Read進行合並。
1） 獲得單倍型 "step1: deal with mutfile and get haplotypes"：varben.deal_mut.checkMutInput.get_haplotypes
① 判悔轎斷突變類型 varben.deal_mut.checkMutInput.check_haplotype -> varben.deal_mut.checkMutInput.check_mut_info
2）獲得序列： "step2: deal haplotypes and get total_chosen_reads, total_chosen_reads_muts"
varben.deal_mut.dealHaplotype.deal_haplotype_multi
1）讀碧族肆取結構變化(SV)模擬的類型分類：
類型1：Read2與斷點A重合或者Read1與斷點B重合
類型2：Read1與斷點A重合或者Read2與斷點B重合
類型3：斷點落在配對Reads中間
類型4: 配對Read完全落在兩個斷點之間
類型5：斷點同時與Read1和Read2重合
2）SV變異模擬：
i）根據BAM文件中指定的VAF，選擇與指定斷點有Overlap的Read；
ii）選擇的Read根據刪除區域分為五種類型，這些Read被編輯以產生 soft-clipper Reads和另一端Read；
iii）編輯後的Reads重新進行比對；
iv）合並重新比對的Read和的原始BAM中的剩餘Read；
0）檢測輸入文件格式 "step0: prepare sv list" ：varben.deal_sv.checkSVInput.check_sv_file
1）獲得插入序列 "step1: get insert size of paired reads"：varben.common.methods.get_insertSize_range
都是 [100, 1000]
2）獲得編輯序列："step2: deal with sv and get total edited reads"
①varben.deal_sv.dealSVType.deal_sv -->
②varben.deal_sv.dealSVType.deal_cnv -->
③varben.deal_sv.getReadsType.pos_type_classify --> ###裡面設置了process=20,
④varben.deal_sv.getReadsType.posType_sub_paired -->
varben.deal_sv.dealSVType.get_write_reads --> ###輸出需要增加或者刪除或者修改的序列ID
3）將需要處理的序列寫入Bam ："step3: get reads by region bed and write bam file"
①varben.deal_sv.getReadsByRegion.get_reads_by_region -->
②varben.common.bamconvert.getRegionReads --> ###生成bam
samtools view raw.bam -b -h
-o cnv_out/tempDir/used_tmp.bam
-U cnv_out/tempDir/exclude_tmp.bam
-L cnv_out/tempDir/consider_region.bed
③varben.common.bamconvert.bamIndex ---> ###創建bam索引
④varben.deal_sv.writeBamByChr.write_sub_bam
4）合並編輯後的read到原始文件中，並重新比對生成新的bam文件："step4: merge edited reads and remap to new bam, consider about the tag, RG, life reads"
varben.deal_sv.mergeEditBam.merge_edit_bam
5）合並bam文件："step5: remapped edit reads and merge"
varben.common.bamconvert.remap --> varben.common.bamconvert.bamToFastq
(edit.bam -> fq1,fq2 -> edit.rmap.bam -> edit.remap.sort.bam -> edit.sort.bam)
Ⅲ DNA/RNA序列比對軟體整理
文章僅是記錄自己的學習使用，有錯誤請指出，我立刻改正
在對比對工具進行比較時，通常將其分為DNA比對工具（DNA-seq）和RNA比對工具（RNA-seq）。它們的區別在於是否會考慮跨外顯子的比對，即：是否會將沒有比對上的reads劈開，對劈開後的兩部分再次比對）。
隨著現在各種seq測序的出現，我們已經不能簡單的根據是比對DNA還是RNA來判斷。比對工具的選擇主要依據reads的比對是否需跨外顯子。（PRO-seq/GRO-seq，它們雖然在建庫時捕獲的RNA，但是它們的比對並不需要考慮跨外顯子。）
常用工具：
DNA-seq：BWA；bowtie&bowtie2
RNA-seq：STAR；HISAT2；Tophat&Tophat2
BWA主要應用二代測序後的大量短小片段與參考基因組之間的定位比對。需要先嫌燃對參考序列建建立索引，BWA也是基於 BWT和 FM-Index 理論來對參考基因組做索引。根據測序方法的不同，有單末端序列(Single-end，SE)比對和雙末端序列(Pair-end，PE)比對。
bowtie出現在測序行業還不成熟的時候，序列長度普遍在50bp以下，bowtie的只滿足長度在50bp以下的reads的比對。官方稱其可以把短的DNA序列（35bp）快速的比對到人類基因組上。
Bowtie2 是一款經典的短讀長序列（ 50-100 bp，最多可到1000 bp ) 比對軟體，節約內存且靈活與成熟的短序列比對軟體，比較適合下一代測序技術。支持單端測序（unpaired) 和雙端測序的比對。支持全局比對（end-to-end align ) 和 局部比對( local align )。其通常使用全文分索引(FM-index)以及Burrows-Wheeler 變換(BWT)索引基因組使得比對非常快速且內存高效，但是這種方法不適合於找到較長的、帶缺口的序列比對
結論：bowtie和bowtie2，是兩個不同類型的比對工具，bowtie2並非是bowtie的升級。尺有所長寸有所短，bowtie適合長度在50b長度以內的reads比對，而bowtie2適合50-100b，甚至更長的reads比對。但是這兩個都屬DNA-seq比對工具
RNA-Seq測序的特性，天然的會有一部分數據延伸到內含子區，這部分跨越外顯子和內含子的reads就稱為『junction reads』，所以RNA-Seq比對軟體需要針對此進行優化。
（ junction：轉錄組reads比對不同於基因組reads比對(如ChIP-seq、WES等)的地方在於，比對的reads可能來源於2個被內含子隔開的外顯子區域，導致reads一端比對在第一個外顯子的後面部分，另一端比對在第二個外顯子的前面部分，芹衡虛即跨剪切位點，從而形成exon-exon junction (剪接點)。這些reads又稱為junction reads，對轉錄本的拼接、鑒定和差異分析具有重要的意義。）
（soft-clip事件: 即reads末端存在低質量鹼基或接頭導致比對不上的， STAR會自動嘗試截去未比對部分，只保留比對上的部分。）
STAR是ENCODE皇家御用的RNA-seq比對工具，ENCODE計劃(ENCyclopedia Of DNA Elements)又稱人類基因組DNA元件網路全書計劃，是2003年在人類基因組計劃完成之後緊接著的又一個大型國際科研項目。
Tophat2的原作者們也不知道是出於什麼考慮，不再更新Tophat2，轉而開發了一個新的比對工具HISAT2，更是推薦人們使用HISAT2，聲稱其速度更快，內存佔用率更小，准確率更高。
此外，HISAT2不僅支持RNA-seq的比對還支持DNA-seq比對，唯一需要做的就是加上一個參數--no-spliced-alignment。但是就目前攔散來看，大部分人都是使用HISAT2做RNA-seq，沒人使用它做DNA-seq
Tophat/Tophat2工具本身不能進行比對，它是通過調用bowtie/bowtie2進行比對的。劃重點，bowtie2不是bowtie的升級版，但是Tophat2是Tophat2的升級版。因此Tophat只可以調用bowtie，而Tophat2不僅可以調用bowtie2（默認）還可以更改設置調用bowtie。
Tophat/Tophat2調用bowtie/bowtie2後，會首先使用bowtie/bowtie2對序列進行比對，對於那些沒有比對上的，會考慮其跨外顯子的可能性，將reads劈開重新比對。
全長轉錄組（Full-length transcriptome）是基於PacBio和Nanopore三代測序平台，無需打斷拼接，直接獲得包含5』UTR、3』UTR、polyA尾的mRNA全長序列及完整結構信息，從而准確分析有參考基因組物種可變剪接及融合基因等結構信息，克服無參考基因組物種轉錄本拼接較短、信息不完整的難題。同時還可以藉助二代測序數據，進行轉錄本特異性表達分析，獲得更加全面的注釋信息。
傳統的使用比較多的長讀長比對軟體是GMAP，05年發表公布，最開始是用來比對低通量的est序列的，後來也有進一步升級為GSNAP支持高通量的二代測序。PacBio測序技術出現後，常用於Iso-seq轉錄本的鑒定，目前仍是相關研究引用量最高的比對軟體，該軟體也一直在持續更新升級。其可以將轉錄本序列與參考基因組序列比對，輸出gff文件，比對速度稍慢。
Minimap2是生信大牛李恆18年用C語言開發的可以用於三代數據（subreads、iso-seq）比對的長序列比對軟體，與傳統的三代比對工具GMAP相比，其速度有非常顯著的提升，當然同時消耗的內存也比較大。使用方法也比較簡單，近幾年引用次數增長的也很迅速，所以大家可以試試用minimap2進行Iso-seq的比對。
Ⅳ 如何利用python進行高通量測序數據的分析 知乎
建議讀兩本書：
1、集體智慧編程 (豆瓣)
因為Python是一門不需要花太多精力（甚至可以說很少），就可以基本掌握的一門語言，所以推薦這本書。題主提到以後想學機器學習，這是一本非常好的入門書，書中的例子源碼都是Python實現的，並且能幫你迅速熟悉Python相關的各種計算庫。
2、統計學習方法 (豆瓣)
考慮到題主要學得踏實，這本書深入淺出地講了和機器學習有關的一切數學基礎知識，一整本的干貨，沒有廢話，非常值得一讀。題主數學專業的話，讀起來應該會比我更順暢。
Ⅳ 以下哪些工具可以獲得序列標識圖
序列標識圖是指將蛋白質序列上的結構和功能信息映射到一個二維圖上，以便更好地理解和分析蛋白質序列的結構和功能。以下是一些可以獲得序列標識圖的工具:
1. InterProScan: InterProScan是一種基於序列的功能注釋工具，可以預測蛋白質序譽察列的結構和功能，並輸出序列標識圖。它使用多種演算法和資料庫，包括Pfam、PRINTS、ProSite、SUPERFAMILY等，以提高預測准確性。
2. Protter: Protter是一種在線的序列標識圖生成工具，它可以將FASTA格式的蛋白質序列映射到一個可視化的二維圖上，顯示出蛋白喊飢質的各種結構和功能信息，包括域、保守區域、跨膜區域、信號肽等。
3. ESPript: ESPript是一種在線的序列標識圖生成工具，它可以將蛋白質序列映射到一個二維圖上，以顯示出蛋白質的二級結構、保守區域、氨基酸殘基間的相互作用慶滲茄等信息。它還支持多種輸出格式和樣式，用戶可以根據需要進行調整。
4. WebLogo: WebLogo是一種在線的序列標識圖生成工具，它可以將多個蛋白質序列的保守區域映射到一個二維圖上，以顯示出每個位置上的氨基酸分布情況，並計算每種氨基酸出現的頻率。它還支持多種輸出格式和樣式，用戶可以根據需要進行調整。
總之，以上工具都可以用於生成蛋白質序列標識圖，用戶可以根據自己的需要選擇合適的工具，並根據輸出結果進行分析和解釋。
Ⅵ ATAC-seq專題---生信分析流程
ATAC-seq信息分析流程主要分為以下幾個部分：數據質控、序列比對、峰檢測、motif分析、峰注釋、富集分析，下面將對各部分內容進行展開講解。
下機數據經過過濾去除接頭含量過高或低質量的reads，得到clean reads用於後續分析。常見的trim軟體有Trimmomatic、Skewer、fastp等。fastp是一款比較新的軟體，使用時可以用--adapter_sequence/--adapter_sequence_r2參數傳入接頭序列，也可以不填這兩個參數，軟體會自動識別接頭並進行剪切。如：
fastp \
--in1 A1_1.fq.gz \ # read1原始fq文件
--out1 A1_clean_1.fq.gz \ # read1過濾後輸出的fq文件
--in2 A1_2.fq.gz \ # read2原始fq文件
--out2 A1_clean_2.fq.gz \ # read2過濾後輸出的fq文件
--cut_tail \ #從3』端向5』端滑窗，如果窗口內鹼基的平均質量值小於設定閾值，則剪切
--cut_tail_window_size=1 \ #窗口大小
--cut_tail_mean_quality=30 \ #cut_tail參數對應的平均質量閾值
--average_qual=30 \ #如果一條read的鹼基平均質量值小於該值即會被舍棄
--length_required=20 \ #經過剪切後的reads長度如果小於該值會被舍棄
fastp軟體的詳細使用方法可參考：https://github.com/OpenGene/fastp。fastp軟體對於trim結果會生成網頁版的報告，可參考官網示例http://opengene.org/fastp/fastp.html和http://opengene.org/fastp/fastp.json，也可以用FastQC軟體對trim前後的數據質量進行評估，FastQC軟體會對單端的數據給出結果，如果是PE測序需要分別運行兩次來評估read1和read2的數據質量。
如：
fastqc A1_1.fq.gz
fastqc A1_2.fq.gz
FastQC會對reads從鹼基質量、接頭含量、N含量、高重復序列等多個方面對reads質量進行評估，生成詳細的網頁版報告，可參考官網示例：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/good_sequence_short_fastqc.html
經過trim得到的reads可以使用BWA、bowtie2等軟體進行比對。首先需要確定參考基因組fa文件，對fa文件建立索引。不同的軟體有各自建立索引的命令，BWA軟體可以參考如下方式建立索引：
bwa index genome.fa
建立好索引後即可開始比對，ATAC-seq推薦使用mem演算法，輸出文件經samtools排序輸出bam：
bwa mem genome.fa A1_clean_1.fq.gz A1_clean_2.fq.gz
| samtools sort -O bam -T A1 > A1.bam
值得注意的是，在實驗過程中質體並不能完全去除，因此會有部分reads比對到質體序列上，需要去除比對到質體上的序列，去除質體序列可以通過samtools提取，具體方法如下：首先將不含質體的染色體名稱寫到一個chrlist文件中，一條染色體的名稱寫成一行，然後執行如下命令即可得到去除質體的bam
samtools view -b A1.bam $chrlist > A1.del_MT_PT.bam
用於後續分析的reads需要時唯一比對且去重復的，bwa比對結果可以通過MAPQ值來提取唯一比對reads，可以用picard、sambamba等軟體去除p，最終得到唯一比對且去重復的bam文件。
比對後得到的bam文件可以轉化為bigWig（bw）格式，通過可視化軟體進行展示。deeptools軟體可以實現bw格式轉化和可視化展示。首先需要在linux環境中安裝deeptools軟體，可以用以下命令實現bam向bw格式的轉換：
bamCoverage -b A1.bam -o A1.bw
此外，可以使用deeptools軟體展示reads在特定區域的分布，如：
computeMatrix reference-point \ # reference-pioint表示計算一個參照點附近的reads分布，與之相對的是scale-regions，計算一個區域附近的reads分布
--referencePoint TSS \#以輸入的bed文件的起始位置作為參照點
-S A1.bw \ #可以是一個或多個bw文件
-R gene.bed \ #基因組位置文件
-b 3000 \ #計算邊界為參考點上游3000bp
-a 3000 \ #計算邊界為參考點下游3000bp，與-b合起來就是繪制參考點上下游3000bp以內的reads分布
-o A1.matrix.mat.gz \ #輸出作圖數據名稱
#圖形繪制
plotHeatmap \
-m new_A1.matrix.mat.gz \ #上一步生成的作圖數據
-out A1.pdf \ # 輸出圖片名稱
繪圖結果展示：
MACS2能夠檢測DNA片斷的富集區域，是ATAC-seq數據call peak的主流軟體。峰檢出的原理如下：首先將所有的reads都向3'方向延伸插入片段長度，然後將基因組進行滑窗，計算該窗口的dynamic λ，λ的計算公式為：λlocal = λBG（λBG是指背景區域上的reads數目），然後利用泊松分布模型的公式計算該窗口的顯著性P值，最後對每一個窗口的顯著性P值進行FDR校正。默認校正後的P值（即qvalue）小於或者等於0.05的區域為peak區域。需要現在linux環境中安裝macs2軟體，然後執行以下命令：
macs2 callpeak \
-t A1.uni.dep.bam \ #bam文件
-n A1 \ # 輸出文件前綴名
--shift -100 \ #extsize的一半乘以-1
--extsize 200 \ #一般是核小體大小
--call-summits #檢測峰頂信息
註：以上參數參考文獻（Jie Wang，et.al.2018.「ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.」Nature Communications）
ATAC分析得到的peak是染色質上的開放區域，這些染色質開放區域常常預示著轉錄因子的結合，因此對peak區域進行motif分析很有意義。常見的motif分析軟體有homer和MEME。以homer軟體為例，首先在linux環境中安裝homer，然後用以下命令進行motif分析：
findMotifsGenome.pl \
A1_peaks.bed \ #用於進行motif分析的bed文件
genome.fa \ #參考基因組fa文件
A1 \ #輸出文件前綴
-size given \ #使用給定的bed區域位置進行分析，如果填-size -100,50則是用給定bed中間位置的上游100bp到下游50bp的區域進行分析
homer分析motif的原理及結果參見：http://homer.ucsd.e/homer/motif/index.html
根據motif與已知轉錄因子的富集情況可以繪制氣泡圖，從而可以看到樣本與已知轉錄因子的富集顯著性。
差異peak代表著比較組合染色質開放性有差異的位點，ChIP-seq和ATAC-seq都可以用DiffBind進行差異分析。DiffBind通過可以通過bam文件和peak的bed文件計算出peak區域標准化的readcount，可以選擇edgeR、DESeq2等模型進行差異分析。
在科研分析中我們往往需要將peak區域與基因聯系起來，也就是通過對peak進行注釋找到peak相關基因。常見的peak注釋軟體有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker為例，需要在R中安裝ChIPseeker包和GenomicFeatures包，然後就可以進行分析了。
library(ChIPseeker)
library(GenomicFeatures)
txdb<- makeTxDbFromGFF(『gene.gtf』)#生成txdb對象，如果研究物種沒有已知的TxDb,可以用GenomicFeatures中的函數生成
peakfile <-readPeakFile(『A1_peaks.narrowPeak』)#導入需要注釋的peak文件
peakAnno <- annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)
# 用peak文件和txdb進行peak注釋，這里可以通過tssRegion定義TSS區域的區間
對於peak注釋的結果，也可以進行可視化展示，如：
p <- plotAnnoPie(peakAnno)
通過注釋得到的peak相關基因可以使用goseq、topGO等R包進行GO富集分析，用kobas進行kegg富集分析，也可以使用DAVID在線工具來完成富集分析。可以通過挑選感興趣的GO term或pathway進一步篩選候選基因。