準備工作：
- 準備數(shù)據(jù)
- 參考基因組：Ler-1.allpaths_lg.final.assembly.fasta
- HiC數(shù)據(jù)：data_1.fastq.gz data_2.fastq.gz
- 安裝所需軟件并軟連接到~/.local下。（添加環(huán)境變量）
- bwa
- samtools
其他支持性軟件： - bedtools
- lachesis
- R
工作流程
建立參考基因組的bwa索引
bwa的比對沒有bowtie2那么嚴格。
mkdir ref && cd ref
bwa index Ler-1.allpaths_lg.final.assembly.fasta
數(shù)據(jù)比對
cd ..
mkdir 02.bwa && cd 02.bwa
bwa mem ../ref/Ler-1.allpaths_lg.final.assembly.fasta -t 10 ../01.fq/data_1.fastq.gz ../01.fq/data_2.fastq.gz > ninanjie.sam
數(shù)據(jù)過濾
- 過濾掉比對時大于2000表示分段匹配結果的sub-alignment。
perl /data/software/3dgenome/pip/LACHESIS/PreprocessSAMs-rmsubalig.pl
ninanjie.sam ninanjie.II.sam
- 過濾距離酶切位點500bp以外的reads，并選取唯一比對的reads。
這一步需要用到PreprocessSAMs.pl。我們需要用vim打開這個腳本修改一些內(nèi)容以適應當前所需要處理的物種。
vim PreprocessSAMs.pl
測試物種是擬南芥，HiC實驗中酶切使用的是HindⅢ，對應的酶切位點序列是AAGCTT，因此需要修改$RE_site = 'AAGCTT'（一般現(xiàn)在用四堿基酶比較多，因為四堿基酶的酶切位點44比六堿基酶的46更多，分辨率更高。）
這里還需要指定Lachesis內(nèi)部的一個腳本
make_bed_around_RE_site.pl的位置還有bedtools和samtools的安裝位置。另外兩個$mem和$picard_head就注釋著不用管。接下來修改
PreprocessSAMs.sh文件
vim PreprocessSAMs.sh
這里需要將
DIR=設置成上一步PreprocessSAMs.pl所在的文件夾，把ASMs=設置成前面生成的ninanjie.II.sam,ASSEMBLY=設置成參考基因組所在的位置。運行
PreprocessSAMs.sh腳本
sh PreprocessSAMs.sh
cd ..
到這里為止前期的數(shù)據(jù)準備階段就完成了，接下來就是激動人心的Lachesis組裝階段了~
mkdir 03.lachesis && cd 03.lachesis
mkdir bam_file
cd bam_file
ln -s /path/to/samplex.II.REduced.paired_only.bam
ln -s /path/to/sampley.II.REduced.paired_only.bam
cd ..
復制一個test_case.ini文件到當前目錄
cp /path/to/test_case.ini
這里需要根據(jù)物種的不同修改很多的參數(shù)。
SPECIES可以不修改不影響結果。OUTPUT_DIR設置成out/90_2_3_120_10，這里后面的這串數(shù)字是下面要設置的各項參數(shù)，建議先把下面設置好之后再回來設置這一項。
DRAFT_ASSEMBLY_FASTA = 修改為參考基因組的路徑SAM_DIR 后面改為比對文件所在的路徑（即上一步的bam_file）SAM_FILES 后面加過濾得到的bam文件的名字（擬南芥的有兩組數(shù)據(jù)samplex.II.REduced.paired_only.bam sampley.II.REduced.paired_only.bam ）RE_SITE_SEQ 后面接酶切序列（擬南芥是AAGCTT）
USE_REFERENCE 0代表不使用標準基因組進行評估1表示用標準基因組進行評估REF_ASSEMBLY_FASTA 后面是標準參考基因組（即該物種最權威的基因組版本）BLAST_FILE_HEAD后面加標準基因組和參考基因組的blastn比對結果（我們自己的參考序列版本與權威版本TAIR10進行BLASTN比對的結果）
這里基本不用修改，如需修改請按照注釋內(nèi)容進行修改。
CLUSTER_N =物種所對應的染色體數(shù)目（擬南芥=5）CLUSTER_MIN_RE_SITES（聚類分析中contig/scaffold序列中的最小酶切位點數(shù)，建議設為90）CLUSTER_MAX_LINK_DENSITY（聚類分析中contig/scaffold序列中的最大Link深度，建議設為2）
CLUSTER_NONINFORMATIVE_RATIO（聚類回插中contig/scaffold序列與目標cluster互作數(shù)同其他cluster互作數(shù)比例，建議設為3。即待回插的序列與目標cluster的互作強度大于其他序列的3倍才予以回插）ORDER_MIN_N_RES_IN_TRUNK（ 排序分析中TRUNK內(nèi)contig/scaffold序列中的最小酶切位點數(shù)，建議設為120。前面設置90的是單條contig/scaffold序列中的最小酶切位點數(shù)，這里是由contig/scaffold連接成TRUNK之后的最小酶切位點數(shù)）
以上這些參數(shù)都是需要根據(jù)不同的物種進行多次調節(jié)嘗試的。
運行l(wèi)achesis
mkdir ~/bin
mkdir ~/public_html
#新建的這兩個目錄是作為組裝基因組和標準基因組比對圖片的存放地址
export PATH=/path/to/R/R-3.4.1/bin/:$PATH
export PATH=/path/to/shendurelab-LACHESIS-c23474f/:$PATH
ulimit -s unlimited
#對上面這一項不熟悉的可以參考→https://zhidao.baidu.com/question/753187924640549364.html
Lachesis test_case.ini
新建的兩個目錄是作為組裝基因組和標準基因組比對圖片的存放地址。
運行需要好一會兒，建議網(wǎng)絡不穩(wěn)定的小伙伴最后一步掛nohup或者開screen。
結果展示
感覺圖還是挺好看的~放到文章里也完全ok。
2018-10-1 updates:
培訓的時候寫的部分記錄，果然有些東西就跟晚上的夢一樣，一定要馬上記錄下來否則隨著時間就漸漸消散得無影無蹤了。
HiC培訓問答拾遺
Q: 如何區(qū)分compartment？
A: 通過看GC含量。compartment A和compartment B會有明顯的GC含量的區(qū)別。從各種圖上都是無法直接得出直觀的結果的，需要依據(jù)矩陣計算。
Q: loop和TAD是什么關系？
A: 他們只是不同的切入點而已。
Q: 有哪些因素會影響HiC的分辨率
A: 主要有三個：測序深度，內(nèi)切酶和bin的大小。
測序深度越深，HiC分辨率越高（當然到后期提升是會逐漸趨于平緩的）。
四堿基內(nèi)切酶分辨率會比六堿基內(nèi)切酶分辨率更高。從概率上講四堿基酶每44個堿基就會出現(xiàn)一個剪切位點，而六堿基酶則每46個堿基才會出現(xiàn)一個酶切位點。
bin相當于畫圖時候的像素點，bin的size選擇得越小，數(shù)量越多，分辨率越高。
目前普遍能接受10kb左右的分辨率。
HiC的理論基礎
- 染色體在核內(nèi)是存在染色體疆域（Chromatin Territory）的。即染色體與染色體之間不是隨意混合交纏在一起的，而是涇渭分明有自己明顯的區(qū)域與界限的。實驗支持：用射線破壞細胞核的某個點，受到損傷的總是部分染色體，而不是所有的染色體都受到損傷。
- 順式作用高于反式作用。順式作用指同一染色體內(nèi)的相互作用，反式作用指染色體之間的相互作用。因為染色體疆域的存在，順式作用會遠遠高于反式作用。
- 互作強度隨線性距離增加而減弱。
