实验目标： （1）首先，正常的，将数据重命名（使用rename.py），建立索引。 （2）比对富集。 （3）试图绘制，行为样本，列为不同染色质标记的热图。 现在开始执行。
awk -v FS="\t" '{print $1,"\t"$3,"\t"$2,"\t"$5,"\t"$6,"\t"$7,"\t"$4,"\t"$1."_sort_peaks.narrowPeak.bed"}' human_hm_full_QC.txt >humdata.txt mv human_hmdata.txt hm_human_information_data.txt python rename.py -m ./hm_human_information_data.txt -i ./human_hm/ -o ./named -n cistrome_id_to_name_map.txt
这个到目前而言是顺利的，得到的结果文件：
接下来，同样的进行排序，这个是建立索引的必要的步骤。 进行排序的时候，出现了这样的问题。
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "./named/[A-J]*" ./named_sort/ 30
>xargs: unmatched single quote; by default quotes are special to xargs unless you use the -0 option H3K27ac_LIS2_Embryonic_Stem_Cell_Embryo_.1.bed H3K27ac_LM7_Osteosarcoma_cell_None_.0.bed H3K27ac_LIS2_Embryonic_Stem_Cell_Embryo_.0.bed ls: write error: Broken pipe
不知道这个报错对最后的结果是否有影响。【会影响】
受文章的启发：https://www.ducea.com/2007/11/22/xargs-unmatched-single-quote/#:~:text=After a while the command failed with the,name and the quote used to quote names.
终于找到了错误的原因：
#!/bin/sh
set -euo pipefail
if [ "$#" -lt "2" ]; then
echo "usage: `basename $0` <input path> <output dir> <threads>"
exit 0
fi
export IN_PATH=$1
export OUT_PATH=$2
export THREADS=1
if [ "$#" -eq "3" ]; then
export THREADS=$3
fi
SORT()
{
set -euo pipefail
IN=$1
if [ ${IN: -4} == ".bed" ]; then
BASE=`basename $IN`
echo $BASE
cat $IN | LC_ALL=C sort --buffer-size 2G -k1,1 -k2,2n -k3,3n | bgzip -c > $OUT_PATH/$BASE.gz
elif [ ${IN: -7} == ".bed.gz" ]; then
BASE=`basename $IN`
echo $BASE
gzip -d -c $IN | LC_ALL=C sort --buffer-size 2G -k1,1 -k2,2n -k3,3n | bgzip -c > $OUT_PATH/$BASE
else
echo "SKIPPING $IN"
fi
}
export -f SORT
ls $IN_PATH | xargs -I{} -t -P $THREADS bash -c "SORT {}" #这里就是错误的代码行
而我对此是这样修改的(加了一个参数-Q)：
ls $IN_PATH -Q | xargs -I{} -t -P $THREADS bash -c "SORT {}"
在ls中，-Q的解释是这样的:
-Q, --quote-name enclose entry names in double quotes
--quoting-style=WORD use quoting style WORD for entry names:
literal, locale, shell, shell-always, c, escape
就是说，我们正式的规定，ls的输出的结果（要传入的参数），必须要用双引号括起来(下面的代码行展示的就是-Q的功能，是我瞎猫碰死耗子碰出来的)。因为下游xargs在进行分割的时候，对于文本字符和非文本字符要求很高，就很容易混淆，发生错误。说来也是好让人生气。困惑了我好久的问题，解决的方法居然如此的简单。
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort -Q |head
"CENPA_HeLa_Epithelium_Cervix_.0.bed.gz"
"CENPA_HeLa_Epithelium_Cervix_.1.bed.gz"
"CENPA_HeLa_Epithelium_Cervix_.2.bed.gz"
"CENPA_HeLa_Epithelium_Cervix_.3.bed.gz"
"CENPA_HT1080_fibrosarcoma_None_.0.bed.gz"
"CENPA_HT1080_fibrosarcoma_None_.1.bed.gz"
"CENPA_HT1080_fibrosarcoma_None_.2.bed.gz"
"CENPA_HuRef_Lymphoblast_None_.0.bed.gz"
"CENPA_HuRef_Lymphoblastoid_None_.0.bed.gz"
"CENPA_HuRef_Lymphoblastoid_None_.1.bed.gz"
然后接下来就是建立索引，这一块是很容易报错的。
面对这个问题，突然想“荡开一笔”，不那么想一口吞一个胖子，想看看分开来，是什么意思（有时候遇到难题的时候，往往是这种瞬间的奇思妙想拯救了我）。
我好像看到了一些希望。这个好像比我想象的要更加的复杂。
我们从中选择一些比较重要的，用于接下来的分析。
awk -v FS="\t" '{print $7}' hm_human_information_data.txt |sort |uniq
H3K27ac ：In enhancers, its presence is strongly anticorrelated with that of the activating histone mark H3k27ac
2539 H3K27ac
H3K27me3：异染色质区域，下调附近基因。
1440 H3K27me3
H3K36me3：基因转录区域。
608 H3K36me3
H3K4me1：增强子。
1210 H3K4me1
H3K4me3：转录激活，启动子。
2313 H3K4me3
H3K9me3：异染色质区域。
590 H3K9me3
在做之前，我们需要验证一下，这些文件是不是都全了。
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort/ -Q |grep 'H3K27ac' |wc -l
2529 #缺10
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort/ -Q |grep 'H3K27me3' |wc -l
157 #缺的最多
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort/ -Q |grep 'H3K36me3' |wc -l
77 #缺的也多
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort/ -Q |grep 'H3K4me1' |wc -l
1207
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort/ -Q |grep 'H3K4me3' |wc -l
2306
(base) [xxzhang@cu08 human_histone_mark]$ ls ./named_sort/ -Q |grep 'H3K9me3' |wc -l
94 #缺的也多
如果有可能，这些缺的是不是那些为0的文件。
重新来。
sed 's/ \t/\t/g' hm_human_information_data.txt >hm_human_information_data_v2.txt
后来经过一系列的探索，发现问题的原因不是“重命名”的文件存在特殊的字符，而原始文件的命名格式，和我们预期的不一样，所以会有一些不符合我们预期规范的成为漏网之鱼。
这个事情告诉我们，很多事情不能想当然，随时检验是多么的重要。
所以，接下来我的任务就是批量的规范化命名目录下的文件。
这里，由于我很弱的shell基础，这里也探索了好久，网上并没有现成的，需要根据一些混乱的代码的提示，自己想办法在此基础上修改。
这里学到的一个很重要的概念是【变量替换】。
以下是我运行成功的shell指令。
#! /bin/bash
for var in `ls *.bed`
do
mv $var ${var%%_*.bed}.bed
done
这个代码的关键在第三行，%%*的意思是从匹配的变量的右边开始删除，知道最后一个""符号。
匹配到的变量譬如是“16790_b_sort.narrowPeak”或“1216_Hsdss_sort”……不管多么的奇奇怪怪，然后删除到从右边开始最后一个_，于是，只剩下数字部分，然后新的文件名就是数字加上“.bed”。
在hm_human_information_data_v2.txt的基础上修改文件，变成hm_human_information_data_v3.txt。
(base) [xxzhang@cu08 human_hm]$ bash rename.sh
(base) [xxzhang@cu08 human_hm]$ mv rename.sh ..
(base) [xxzhang@cu08 human_hm]$ cd ..
(base) [xxzhang@cu08 human_histone_mark]$ mv named named_pre
(base) [xxzhang@cu08 human_histone_mark]$ mkdir named
(base) [xxzhang@cu08 human_histone_mark]$ python rename.py -m ./hm_human_information_data_v3.txt -i ./human_hm/ -o ./named -n cistrome_id_to_name_map_v3.txt
(base) [xxzhang@cu08 human_histone_mark]$ cd named
(base) [xxzhang@cu08 named]$ls |wc -l
11060
(base) [xxzhang@cu08 named]$./sort_bed "./named/*" ./named_sort/ 30
#然后成功了。
接下来，对其进行排序，文件数量基本上一致(11016/11060，去掉一些空文件)。
(base) [xxzhang@cu08 named_sort]$ ls |wc -l
11016
(base) [xxzhang@cu08 named_sort]$ cd ..
(base) [xxzhang@cu08 human_histone_mark]$ cd named
(base) [xxzhang@cu08 named]$ ls |wc -l
11060
然后接下来，我们再同一个我们感兴趣的那几个marker的文件的数量。
(base) [xxzhang@cu08 human_histone_mark]$ cd named_sort
(base) [xxzhang@cu08 named_sort]$ ls |grep 'H3K27ac' |wc -l
2528/2539
(base) [xxzhang@cu08 human_histone_mark]$ ls |grep 'H3K27me3' |wc -l
1433/1440
(base) [xxzhang@cu08 human_histone_mark]$ ls |grep 'H3K36me3' |wc -l
611/608
(base) [xxzhang@cu08 human_histone_mark]$ ls |grep 'H3K4me1' |wc -l
1206/1210
(base) [xxzhang@cu08 human_histone_mark]$ ls |grep 'H3K4me3' |wc -l
2302/2313
(base) [xxzhang@cu08 human_histone_mark]$ ls |grep 'H3K9me3' |wc -l
588/590
- H3K27ac [暂时不行]
(base) [xxzhang@cu08 human_histone_mark]$ time giggle index -i "./named_sort/H3K27ac*" -o ./named_sort_H3K27ac_b -s -f
Indexed 34094652 intervals.
real 5m18.042s
user 3m5.982s
sys 0m14.399s
这里成功了，似乎有了一点点的希望。
接下来按照这个思路，继续去构建其他的染色质标记的索引。
- H3K4me1 [暂时不行]
(base) [xxzhang@cu08 human_histone_mark]$ time giggle index -i "./named_sort/H3K4me1*" -o ./named_sort_H3K4me1_b -s -f
Could not open file './named_sort/H3K4me1_None_T_Lymphocyte_Blood.2.bed.gz'
giggle: Could not open ./named_sort/H3K4me1_None_T_Lymphocyte_Blood.2.bed.gz.
real 0m3.893s
user 0m0.033s
sys 0m0.061s
而原先的文件有差不多10000多个，也就是说在sort的时候，有很多文件被剔除了。这就成了一个新的问题。
所以，这里回到前面，解决前面的问题。【已解决】
但是又出现了can not open的问题，猜测问题是多处理的文件太多太大，内存消耗超过上限。
- H3K4me3 [暂时不行]
- H3K27me3 [暂时不行]
- H3K36me3
(base) [xxzhang@cu08 human_histone_mark]$ time giggle index -i "./named_sort/H3K36me3*" -o ./named_sort_H3K36me3_b -s -f
Indexed 8842032 intervals.
real 1m26.186s
user 0m41.846s
sys 0m4.480s
- H3K9me3
(base) [xxzhang@cu08 human_histone_mark]$ time giggle index -i "./named_sort/H3K9me3*" -o ./named_sort_H3K9me3_b -s -f
Indexed 10522610 intervals.
real 1m41.868s
user 0m49.491s
sys 0m4.980s
所以现在把那个repeat的文件先拿出来，然后用其与我们已经建好的这两个索引去比较。
(base) [xxzhang@cu08 human_histone_mark]$ giggle search -i ./named_sort_H3K36me3_b/ -q Hs_repeat.bed.gz -s >Hs_repeat.bed.gz.giggle.H3K36me3.result
(base) [xxzhang@cu08 human_histone_mark]$ awk '$8>0' Hs_repeat.bed.gz.giggle.H3K36me3.result >repeat_positive.H3K36me3.result
(base) [xxzhang@cu08 human_histone_mark]$ giggle search -i ./named_sort_H3K9me3_b/ -q Hs_repeat.bed.gz -s >Hs_repeat.bed.gz.giggle.H3K9me3.result
(base) [xxzhang@cu08 human_histone_mark]$ awk '$8>0' Hs_repeat.bed.gz.giggle.H3K9me3.result >repeat_positive.H3K9me3.result
可以了。
我今天先到这里，还算是有收获的。
