【实验记录】2022年7月20日-7月22日
作者：互联网
今日目标：（1）阅读文献：GIGGLE；（2）将其应用到cistromDB转录因子的那套数据中；
一、文献整理
1、GIGGLE用以解决什么问题？
结合已有的表观组的数据，如open chromatin,enhancers和transcribed regions,对我们感兴趣的基因组位点进行注释。
2、GIGGLE的优势：
（1）检索的速度快；（2）对结果进行排名，找到最相关的query intervals；
3、上述优势如何实现：
（1）快：使用unified index来确定query和已知的注释文件间的overlap；使用B+树的结构来最小化磁盘空间；
（2）rank：使用Monte Carlo（MC）进行模拟；
4、输出结果：
-log2(P value)和log2(odds ratio)
注：这里的P value和odds ratio是如何得到的？——构建了一个双尾Fisher's 精确性检验的2×2表。
| the number of the intervals that ae in both the query and index file | solely in the query file |
| neither the query file nor the indexed file |
estimated by the difference between the union of the two sets and the quotient pf the mean interval size of both sets and the genome size. |
二、代码尝试
1、giggle的安装
git clone https://github.com/ryanlayer/giggle.gitcd giggle/
make
#接着配置环境变量，即可使用giggle
giggle, v0.6.3
usage: giggle <command> [options]
index Create an index
search Search an index
2、giggle的运行
（1）建库
首先在建库之前，需要对要index的文件进行预处理：
下面是我们要建库的文件，为组蛋白的结合位点。
以下为文件的内容：
chr3 93470333 93470834 peak350929 2966 . 26.36807 306.13461 296.64355 172
chrKI270466.1 894 1208 peak603579 2942 . 42.38830 301.80423 294.21622 174
chr20 50191112 50192112 peak313970 2609 . 19.62702 268.04147 260.96368 793
chr15 20276358 20276645 peak179139 2510 . 28.56448 257.99652 251.06773 145
chrKI270732.1 1289 1839 peak604004 2480 . 23.20901 254.93692 248.03578 393
chr7 58090994 58091552 peak522566 2306 . 26.20435 237.37437 230.60510 414
chrKI270736.1 172268 172555 peak604062 2285 . 25.23166 235.32509 228.57751 144
chr10 125404964 125405256 peak74441 2266 . 37.07228 233.42575 226.69434 152
chrY 56763383 56763655 peak624507 2174 . 41.71397 224.06732 217.41508 136
chrGL000009.2 200751 201160 peak603104 2121 . 33.73744 218.70740 212.10054 142
其实这些文件有的数据真的很少，那么为了简化下一步的步骤，我们首先对数据进行筛选。
a. 删除那些空文件
find . -name "*" -type f -size 0c | xargs -n 1 rm -f
目前用我们的数据出现了一些问题。
giggle index -i /home/xxzhang/data/Epigenome/cistrome/human_factor/human_factor/*.bed -o human_factor_index/ -f
File type not supported '/home/xxzhang/data/Epigenome/cistrome/human_factor/human_factor/1007_sort_peaks.narrowPeak.bed'.
giggle: Could not open /home/xxzhang/data/Epigenome/cistrome/human_factor/human_factor/1007_sort_peaks.narrowPeak.bed.
移除该文件后，问题依旧不得解。
于是，打算先使用作者的示例数据来看。
b. 使用实例数据运行
wget https://github.com/brentp/gargs/releases/download/v0.3.9/gargs_linux
chmod +x gargs_linux #给该文件可执行权限
mv gargs_linux gargs
vim ~/.bashrc #配置路径
source ~/.bashrc
(base) [xxzhang@mu02 giggle_workplace]$ gargs
gargs 0.3.9
Usage: gargs [--procs PROCS] [--sep SEP] [--nlines NLINES] [--retry RETRY] [--ordered] [--verbose] [--stop-on-error] [--dry-run] [--log LOG] COMMAND
mkdir repeat
url="http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/rmsk.txt.gz"
curl -s $url | gunzip -c | cut -f 6,7,8,11,12,13 > repeat/rmsk.bed
url="http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/simpleRepeat.txt.gz"
curl -s $url | gunzip -c | cut -f 2,3,4,17 > repeat/simpleRepeat.bed
url="http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/microsat.txt.gz"
curl -s $url | gunzip -c | cut -f 2,3,4 > repeat/microsat.bed
url="http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/genomicSuperDups.txt.gz"
curl -s $url | gunzip -c | cut -f 2,3,4,5 > repeat/genomicSuperDups.bed
url="http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/chainSelf.txt.gz"
curl -s $url | gunzip -c | cut -f 3,5,6,7,10,11 > repeat/chainSelf.bed
mkdir repeat_sort
giggle/scripts/sort_bed "repeat/*.bed" repeat_sort 4
giggle index -i "repeat_sort/*gz" -o repeat_sort_b -f -s
giggle search -i repeat_sort_b -r 1:200457776-200457776 -f rmsk,simple -v
chr1 200457488 200457811 L2a LINE L2 repeat_sort/rmsk.bed.gz
现在一行行理解上述代码。
#1、curl 这个命令的意思我理解就是向url发送下载文件的请求；该命令后面参数若干，我也记不住；#似乎这种方法并不需要直接下载文件到本地，而是直接处理。又是极简的一种思路了！
#2、gunzip -c 解压缩文件
#3、cut -f 剪切，且剪切时按照fields定位（这里的field我理解的就是空格）。
#4、sort_bed这行代码，直接把文件夹下的bed文件全部排序并压缩了 只是不知道4 这里是什么意思？（查了一下，是线程的意思）
(base) [xxzhang@mu02 repeat]$ /home/xxzhang/workplace/software/giggle/scripts/sort_bed
usage: sort_bed <input path> <output dir> <threads>
#5、通过与示例代码进行比较，我发现会不会问题的原因在以下方面：
#（a）首先我没有排序，即使用sort_bed文件；
#（b）其次我在建立索引的时候，没有加引号；
下面就是处理后的文件，我们发现基本上都是一些重要位置的坐标信息：前三列分别为染色体号，开始，结束。
注：这个simpleRepeat.bed文件倒是挺有意思的，可以把序列也一起显示出来，或许我也可以这样放到我的那个汇总的gene list文件中。
c. 使用自己的数据在反馈上述结果的基础上，再进行测试。
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "/home/xxzhang/data/Epigenome/cistrome/human_factor/human_factor/*.bed" human_factor_sort 32
giggle index -i "human_factor_sort/*.gz" -o human_factor_index -f -s
Indexed 235577307 intervals.
giggle search -i human_factor_index -r 9:2843571-2843877 -v
得到了bed文件的一些结果，所以上述实验成功。
也可以反推到之前失败的原因就是自己总结的那个。
d. 将人特异性的repeat element的intervals的文件作为输入。
我们还可以把文件作为输入来看。
首先构建输入文件。
cat /home/xxzhang/data/Sequencing/multiz_alignment/insertion_backup/repeat_result.txt |awk '{print $1,$4,$5,$6,$7,$8}'>human_specifc_repeat.bed
对其进行压缩。
gzip -9 human_specifc_repeat.bed
然后将该文件作为输入。
(base) [xxzhang@mu02 giggle_workplace]$ giggle search -i human_factor_index -q human_specifc_repeat.bed.gz -v -o |head
Not a BGZF file 'human_specifc_repeat.bed.gz'
giggle: Error loading query file human_specifc_repeat.bed.gz.
似乎输入文件的压缩格式不对，现在出现报错。
现在换用bgzip进行压缩。
bgzip human_specifc_repeat.bed
这次感觉对了，没有报错。
giggle search -i human_factor_index -q human_specifc_repeat.bed.gz -v -o |head
然后进行search的时候，却发现结果为空。
所以，问题又再次在哪里？是不是输入文件的格式不对？
我们看一下作者在实例中，用到的测试数据。这个数据能够运行起来。
比较一下两个输入文件的格式，感觉没有什么区别。难道是其排序且压缩的缘故嘛？
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "/home/xxzhang/workplace/project/giggle_workplace/*.bed" ./ 32
重新压缩再运行，依旧空文件。
甚至换用重复序列的索引，依旧为空（这个按理说一定会有结果的，因为我的输入文件就是重复序列的索引）。
(base) [xxzhang@mu02 giggle_workplace]$ giggle search -i human_factor_index -q human_specifc_repeat.bed.gz -v -o |head
(base) [xxzhang@mu02 giggle_workplace]$ giggle search -i repeat_sort_b -q human_specifc_repeat.bed.gz -v -o |head
所以，问题出在哪里了呢？
如果换用输入文件中的单个坐标，会有结果吗？看一下是真的没有，还是格式的问题。
(base) [xxzhang@mu02 giggle_workplace]$ giggle search -i human_factor_index -r 13:24136818-24137129 -v
chr13 24136956 24137184 peak42704 57 . 4.37756 7.61568 5.75469 23 human_factor_sort/85217_sort_peaks.narrowPeak.bed.gz
对呀，是可以找到结果的。
所以，问题的关键还在于，输入文件的格式的问题。
那么究竟什么样的输入文件的格式才是对的呢？这里可能需要回到原始的代码中来看。
(base) [xxzhang@mu02 giggle_workplace]$ giggle search
giggle, v0.6.3
usage: giggle search -i <index directory> [options]
options:
-i giggle index directory
-r <regions (CSV)>
-q <query file> #这里并没有具体的注释，输入文件的格式。
-o give reuslts per record in the query file (omits empty results)
-c give counts by indexed file
-s give significance by indexed file (requires query file)
-v give full record results
-f print results for files that match a pattern (regex CSV)
-g genome size for significance testing (default 3095677412)
-l list the files in the index
使用vim的set list指令查看，还是能够发现一些不同的。我发现我自己生生的这个文件中有多余的空格，也许问题就在这个地方。
我现在想办法看如何去除这个空格。
使用sed可以嘛？
sed 's/\s\+/\t/g' test.bed >test2.bed #将其中的多个空格，替换为tab键。
我觉得现在我成功了。但是不知道以此作为输入文件能不能得到结果。
继续对其进行压缩，然后处理。
bgzip -c test2.bed >test2.bed.gz
time giggle search -i human_factor_index -q test2.bed.gz -v -o |head
太棒了，我成功了！debug的过程真的是特别考验自己解决问题的能力的过程，我觉得我在这一方面的训练还是很可以的。遇到一个具体的问题，知道往哪个方向去想，去考虑。
前面的过程，都是自己为了跑通这个pipeline所做的测试，接下来就是正式的带入repeat.bed的注释数据，看看结果。如果这里运行成功的话，就把之前的测试的数据全部整理好，删除（最大化的达到极简的目的）。
(base) [xxzhang@mu02 giggle_workplace]$ sed 's/\s\+/\t/g' human_specifc_repeat.bed >human_specifc_repeat_v2.bed
(base) [xxzhang@mu02 giggle_workplace]$ vim human_specifc_repeat_v2.bed
(base) [xxzhang@mu02 giggle_workplace]$ bgzip -c human_specifc_repeat_v2.bed >human_specifc_repeat_v2.bed.gz
(base) [xxzhang@mu02 giggle_workplace]$ time giggle search -i human_factor_index -q human_specifc_repeat_v2.bed.gz -s > human_specifc_repeat.result
real 0m13.427s
user 0m10.042s
sys 0m0.926s
(base) [xxzhang@mu02 giggle_workplace]$ time giggle search -i human_factor_index -q human_specifc_repeat_v2.bed.gz -v -o > human_specifc_repeat_intervals.result
real 9m0.734s
user 2m4.403s
sys 0m31.394s
现在删除之前测试的时候用到的数据文件。
(base) [xxzhang@mu02 giggle_workplace]$ rm -R repeat_sort
(base) [xxzhang@mu02 giggle_workplace]$ rm -R repeat_sort_b/
(base) [xxzhang@mu02 giggle_workplace]$ rm -R repeat
(base) [xxzhang@mu02 giggle_workplace]$ rm test*
(base) [xxzhang@mu02 giggle_workplace]$ rm GSM1218850_MB135DMMD.peak.*
(base) [xxzhang@mu02 giggle_workplace]$ rm human_specifc_repeat.bed
(base) [xxzhang@mu02 giggle_workplace]$ rm human_specifc_repeat.bed.gz
为了继续简化，将v2的文件命名回原来的。
(base) [xxzhang@mu02 giggle_workplace]$ mv human_specifc_repeat_v2.bed human_specifc_repeat.bed
(base) [xxzhang@mu02 giggle_workplace]$ mv human_specifc_repeat_v2.bed.gz human_specifc_repeat.bed.gz
e. 理解该程序的输出结果。
这个结果中包含（汇总信息）每一个文件中一共有多少peak区域，以及在这些peak区域中，匹配到我们感兴趣的区间的peak的数量。在这个过程中一般不会出现计数重复，因为基因组上位置的坐标是唯一的。可能有的文件中本身reads数很多，包括了了一些contig的位点，可能会影响最终的计算的结果。
这里有两个关键的值：odds ratio和fishers two tail的P值；==>我的意思大概是说这个转录因子和我们感兴趣的这个区间，结合的富集的程度；通过这个信息，我们也许会得到比较有意思的结果。（我自己还蛮期待的）
这个文件的意思是，我们感兴趣的那个位点，都是哪些peak匹配上了，并且分别对应的是哪个转录因子的peak文件，通过对这个文件进行分析，我们可以总结提炼出哪些类型的转录因子会偏好的结合在不同的家族上，这里可以做一些比较有意思的研究。
d. 绘制有生物学意义的图（回到文献，看他们那张图能够说明什么，给我们什么启示）。
作者这里比较的是一个两维的数据，一个是有H3K9me3标记的LTR，一个是没有H3K9me3标记的LTR，主要想知道哪些转录因子特异性的和这两种不同的片段结合（非此即彼的关系）。
但是在咱们的这套数据中，我们只是一维的数据。不过我突然想到，也不是不可以有另一维度的数据，即那些我们认为非人特异性插入的重复序列家族（或者说是全部的重复序列的区间）。我想知道，有哪些转录因子是特异性的和这些人特异性的片段结合的。如果能够得到这一点，这可能就是全新的信息。
好的，我们可以先拿全部的试试看。
#全部的重复序列的片段
(base) [xxzhang@mu02 giggle_workplace]$ cp ./
(base) [xxzhang@mu02 giggle_workplace]$ sed 's/\s\+/\t/g' hg38_stat_SINE_LINE_LTR_SVA_v2.txt >all_repeat.bed
(base) [xxzhang@mu02 giggle_workplace]$ vim all_repeat.bed
(base) [xxzhang@mu02 giggle_workplace]$ bgzip -c all_repeat.bed >all_repeat.bed.gz
(base) [xxzhang@mu02 giggle_workplace]$ time giggle search -i human_factor_index -q all_repeat.bed.gz -s > all_repeat.result
real 2m49.920s
user 2m26.874s
sys 0m10.522s
(base) [xxzhang@mu02 giggle_workplace]$ time giggle search -i human_factor_index -q all_repeat.bed.gz -v -o > all_repeat_intervals.result #这个文件略大，感觉要跑一会儿#
#这里的话，基本上把全部的也跑完了，接下来可能就是提取这两个文件中的odds_ratio的数据，画一个和作者类似的dotplot。
#我觉得我在等代码运行结果的时候，与其摸鱼，不如做一些有价值的事情，比如整理自己，看文献，写日记。
#好的，我现在准备做这样的事情。
22:08 代码还是没有跑完，但是提交到服务器上了。
现在尝试提取并合并两个文件中的odds_ratio的值。
#(base) [xxzhang@mu02 giggle_workplace]$ awk '{print $1,$4}' human_specifc_repeat.result |head >./odds_ratio_plot/human_specific_odds.txt #这里出错了
awk '{print $1,$4}' human_specifc_repeat.result >./odds_ratio_plot/human_specific_odds.txt
#(base) [xxzhang@mu02 giggle_workplace]$ awk '{print $1,$4}' all_repeat.result |head >./odds_ratio_plot/all_odds.txt
awk '{print $1,$4}' all_repeat.result >./odds_ratio_plot/all_odds.txt
(base) [xxzhang@mu02 giggle_workplace]$ cd odds_ratio_plot/
(base) [xxzhang@mu02 odds_ratio_plot]$ sed 's/human_factor_sort\///g' all_odds.txt | sed 's/_sort_peaks.narrowPeak.bed.gz//g' >all_odds_v3.txt #通道真好用，可以去掉一些东西#又是一个极简的思路
(base) [xxzhang@mu02 odds_ratio_plot]$ sed 's/human_factor_sort\///g' human_specific_odds.txt | sed 's/_sort_peaks.narrowPeak.bed.gz//g' >human_specific_odds_v3.txt #极简是最优雅的一条解决方案
(base) [xxzhang@mu02 odds_ratio_plot]$ R
R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> all_data<-read.table("all_odds_v3.txt")
> human_specific<-read.table("human_specific_odds_v3.txt")
> colnames(all_data)<-c("TF","odds_all")
> colnames(human_specific)<-c("TF","odds_Hs")
> merge_data<-merge(all_data,human_specific)
> head(merge_data)
TF odds_all odds_Hs
1 1006 0.47956931 0.13710985
2 1007 0.19325949 0.12020182
3 1009 0.02118476 0.83990674
4 1010 0.60468109 0.13898023
5 1011 0.23594378 0.04056846
6 1012 0.26688225 0.12174564
> write.table(merge_data,"odds.txt",quote=F,row.names=F)
odds文件构建好之后，接下来就想与qc的文件注释合并，看一下每一个文件都对应的是哪些转录因子？
> qc_data<-read.table("human_factor_full_QC.txt",sep="\t",fill=T,header=T,na.strings="",comment.char="#",quote="")
> head(qc_data)
DCid Species GSMID Factor Cell_line Cell_type Tissue_type FastQC
1 1 Homo sapiens GSM448026 BTAF1 HeLa Epithelium Cervix 9
2 2 Homo sapiens GSM448027 GAPDH HeLa Epithelium Cervix 9
3 4 Homo sapiens GSM540707 EGR1 K562 Erythroblast Bone Marrow 9
4 6 Homo sapiens GSM460127 TCF4 LS174T Epithelium Colon 9
5 8 Homo sapiens GSM460125 TCF4 LS174T Epithelium Colon 9
6 9 Homo sapiens GSM460124 TCF4 LS174T Epithelium Colon 9
UniquelyMappedRatio PBC PeaksFoldChangeAbove10 FRiP
1 0.0502 0.880 406 0.094809301
2 0.0274 0.931 248 0.175196234
3 0.0845 0.893 500 0.036462381
4 0.0007 0.819 7 0.021363742
5 0.0583 0.346 22 0.004400244
6 0.1380 0.372 1043 0.050072836
PeaksUnionDHSRatio
1 0.3523726
2 0.2470817
3 0.5229698
4 0.3333333
5 0.5102041
6 0.8650915
> odds<-read.table("odds.txt")
> head(odds)
V1 V2 V3
1 TF odds_all odds_Hs
2 1 2.43309999969923 3.91180734283374
3 2 32.2100357010977 4.55942577377889
4 4 0.304161474441714 0.560265117134011
5 6 2.81268582755203 13.2220390009941
6 8 25.958534422235 10.8848430321593
> odds<-read.table("odds.txt",header=T)
> head(odds)
TF odds_all odds_Hs
1 1 2.4331000 3.9118073
2 2 32.2100357 4.5594258
3 4 0.3041615 0.5602651
4 6 2.8126858 13.2220390
5 8 25.9585344 10.8848430
6 9 0.3630817 0.3967151
> colnames(odds)[1]
[1] "TF"
> colnames(odds)[1]<-"DCid"
> merge_data<-merge(qc_data,odds)
> head(merge_data)
DCid Species GSMID Factor Cell_line Cell_type Tissue_type FastQC
1 1 Homo sapiens GSM448026 BTAF1 HeLa Epithelium Cervix 9
2 2 Homo sapiens GSM448027 GAPDH HeLa Epithelium Cervix 9
3 4 Homo sapiens GSM540707 EGR1 K562 Erythroblast Bone Marrow 9
4 6 Homo sapiens GSM460127 TCF4 LS174T Epithelium Colon 9
5 8 Homo sapiens GSM460125 TCF4 LS174T Epithelium Colon 9
6 9 Homo sapiens GSM460124 TCF4 LS174T Epithelium Colon 9
UniquelyMappedRatio PBC PeaksFoldChangeAbove10 FRiP
1 0.0502 0.880 406 0.094809301
2 0.0274 0.931 248 0.175196234
3 0.0845 0.893 500 0.036462381
4 0.0007 0.819 7 0.021363742
5 0.0583 0.346 22 0.004400244
6 0.1380 0.372 1043 0.050072836
PeaksUnionDHSRatio odds_all odds_Hs
1 0.3523726 2.4331000 3.9118073
2 0.2470817 32.2100357 4.5594258
3 0.5229698 0.3041615 0.5602651
4 0.3333333 2.8126858 13.2220390
5 0.5102041 25.9585344 10.8848430
6 0.8650915 0.3630817 0.3967151
> subData<-merge_data[,c(1,4,6,7,14,15)]
> head(subData)
DCid Factor Cell_type Tissue_type odds_all odds_Hs
1 1 BTAF1 Epithelium Cervix 2.4331000 3.9118073
2 2 GAPDH Epithelium Cervix 32.2100357 4.5594258
3 4 EGR1 Erythroblast Bone Marrow 0.3041615 0.5602651
4 6 TCF4 Epithelium Colon 2.8126858 13.2220390
5 8 TCF4 Epithelium Colon 25.9585344 10.8848430
6 9 TCF4 Epithelium Colon 0.3630817 0.3967151
> write.table(subData,"metadata.txt",quote=F,row.names=F,sep="\t")
接下来这两个数据准备好之后，开始画图（模仿那篇文献中的图）。
(画图的细节明天再看，感觉东西要出来了！然后接下来就继续啃一些细节的问题)
（还有就是继续看文献，整理gene list的数据）
2022.07.22 14:10 在想是怎么回事
(base) [xxzhang@mu02 giggle_workplace]$ grep "73092" human_specifc_repeat.result
human_factor_sort/73092_sort_peaks.narrowPeak.bed.gz 5 4 335.87713655578017 1.481888494486196e-08 0.99999999997797281 1.481888494486196e
(base) [xxzhang@mu02 giggle_workplace]$ grep "73092" all_repeat.result
human_factor_sort/73092_sort_peaks.narrowPeak.bed.gz 5 5 1.5065803663953232e-06 1 1 1 -0
看注释文件上写得是，这边第二列指的是file size，而第三列是overlaps;我原以为指的是有多少片段是位于我们要搜寻的序列中，但是为啥一个是4，一个是5，它的odds_ratio相差那么大呢？因此，肯定是我理解的不太对。
我理解的是对的，但是我的序列中有重复的片段，可能为了得到更加精确的结果，需要重新做一遍。
(base) [xxzhang@mu02 giggle_workplace]$ cat human_specifc_repeat.bed |sort |uniq -c |awk '$1>1'
2 chr10 54757277 54757580 AluYg6 Alu SINE
2 chr1 107922671 107922978 AluYh3 Alu SINE
2 chr11 110974095 110974409 AluYb8 Alu SINE
2 chr11 122510721 122511034 L1PA2 L1 LINE
2 chr11 16466649 16466954 AluYa5 Alu SINE
2 chr11 20650656 20650965 AluY Alu SINE
2 chr1 144924459 144924754 AluYc Alu SINE
2 chr1 156327262 156328873 SVA_D SVA SVA
2 chr1 161959403 161959694 AluYa5 Alu SINE
2 chr1 184309605 184309913 AluYa5 Alu SINE
我不明白，为什么会出现重复的内容呢？#这个会影响我之前的结果吗？
好吧，就在歌声的陪伴下，我重新做一遍。07.22
(base) [xxzhang@mu02 giggle_workplace]$ mkdir new_trial
(base) [xxzhang@mu02 giggle_workplace]$ cat human_specifc_repeat.bed |sort |uniq >./new_trial/human_specific_repeat_uniq.bed
(base) [xxzhang@mu02 giggle_workplace]$ cd new_trial/
(base) [xxzhang@mu02 new_trial]$ vim human_specific_repeat_uniq.bed
(base) [xxzhang@mu02 new_trial]$ bgzip -c human_specific_repeat_uniq.bed >human_specific_repeat_uniq.bed.gz
(base) [xxzhang@mu02 new_trial]$ time giggle search -i ../human_factor_index -q human_specific_repeat_uniq.bed.gz -s > human_specifc_repeat_uniq.result
real 0m32.606s
user 0m9.975s
sys 0m1.091s
出现以下这种情况的主要原因是，存在一对多的这种关系。即在一个chip的peak中，有两段序列，所以，我们能够看到重复的片段。
(base) [xxzhang@mu02 giggle_workplace]$ grep '73092_sort_peaks' all_repeat_intervals.result
chr20 26080324 26080900 peak2 130 . 11.02002 19.28667 13.03425 286 human_factor_sort/73092_sort_peaks.narrowPeak.bed.gz
chr20 26080324 26080900 peak2 130 . 11.02002 19.28667 13.03425 286 human_factor_sort/73092_sort_peaks.narrowPeak.bed.gz
chr20 35975051 35975475 peak3 258 . 15.12442 33.05297 25.81777 219 human_factor_sort/73092_sort_peaks.narrowPeak.bed.gz
chr20 35975051 35975475 peak3 258 . 15.12442 33.05297 25.81777 219 human_factor_sort/73092_sort_peaks.narrowPeak.bed.gz
所以，对于这个peak的重复我们是能够理解的。我想知道，面对这种情况，其计数是什么样的？
后来绘制了这张图，但是这个结果让我再次陷入了迷惑当中。不过这一环节，则基本上告一段落了。
参考文献：[1] Layer RM, Pedersen BS, DiSera T, Marth GT, Gertz J, Quinlan AR. GIGGLE: a search engine for large-scale integrated genome analysis. Nat Methods. 2018 Feb;15(2):123-126.
[2] https://github.com/ryanlayer/giggle
标签：giggle,repeat,20,2022,22,bed,base,xxzhang,human 来源： https://www.cnblogs.com/zjuer/p/16507965.html