0.僅作個人筆記使用
意思就是只考慮我能不能看懂
fithic的安裝、使用這里略去
假定已經(jīng)運行完了fithic
輸入文件長這樣:
六列, 兩部分; 表示顯著互作的兩個區(qū)間
#Start
awk '{if($1!=$4)print}' XX.10k.fithic.chr.bed > XX.10k.fithic.tijian.bed
#處理輸入文件
#得到 染色體體間互作的信息
cut -f1,2,3 XX.10k.fithic.tijian.bed > tmp1
cut -f4,5,6 XX.10k.fithic.tijian.bed > tmp2
cat tmp1 tmp2 | bedtools sort -i - | uniq -c > tmp3
awk '{print $2"\t"$3"\t"$4"\t"$1}' tmp3 > XX.10k.fithic.tijian.density
#X.10k.fithic.tijian.density
#最后的輸出文件, 這個文件后文還會用到
rm tmp1 tmp2 tmp3
#刪除中間文件
輸出文件:
結(jié)果文件長這樣
最后畫一個核型圖
核型圖的繪制參考: https://blog.csdn.net/u013429737/article/details/116121440
最終結(jié)果
2.上述這個圖不能滿足我的需要, 我應(yīng)該是需要meta-plot的樣式:
ref: https://doi.org/10.1038/s41467-021-27091-0
文章的Fig. S16
#Start
#獲取基因組每條染色體的長度
bioawk -c fastx '{print $name,length($seq)}' XX.Chr.fa > XX.length
#拿到每個顯著互作在染色體上的相對位置
#相對位置是這樣定義的:
#該顯著互作所在的位置占整條染色體長度的比例
#比如互作位于染色體的末端, 其位置也就是染色體長度, 相對位置也就是1
python Find_A_in_B.v2.py XX.10k.fithic.tijian.density XX.length | awk '{print $4/$2"\t"$5/$2"\t"$6}' > tmp1
#腳本見教程最后
#此時tmp1的前兩列就是相對位置, 為了更標準化, 我們保留兩位小數(shù)或者三位小數(shù)
#awk '{for(i = 1; i <= NF; i++) {printf("%.2f\t", $i)} {printf("\n")}}' tmp1 > tmp2 #保留兩位小數(shù)
#awk '{for(i = 1; i <= NF; i++) {printf("%.3f\t", $i)} {printf("\n")}}' tmp1 > tmp2 #保留三位小數(shù)
#上述命令行作廢弃甥，因為保留小數(shù)的同時它還做了四舍五入
awk '{printf "%.3f %.3f %.3f\n", int($1*1000)/1000, int($2*1000)/1000, int($3*1000)/1000}' tmp1 > tmp2
#awk中的1000就是保留三位小數(shù)的意思
awk '{print $1"\t"$1+0.001"\t"$3}' tmp2 > tmp3
rm tmp1 tmp2
#獲取bed文件
tmp3的文件內(nèi)容此時是這樣的, 前兩列是相對位置, 第三列
#續(xù)上文
#去重并累加
awk -F "\t" '{sum[$1"\t"$2]+=$3}END{for(c in sum){print c,sum[c]}}' tmp3 | sort -Vk 1 > tmp4
#一定要排序
#awk '{print $2-$1}' tmp4 | sort -u
#如果只輸出一個數(shù)字0.01, 就表明沒問題
awk '{print $3}' tmp4 | tr "\n" "," > XX.meta.tijian.plot.csv
rm tmp3 tmp4
#接下來是R里面的代碼
library(pheatmap)
a = read.csv("XX.meta.tijian.plot.csv", header = F)
#標準化
b = t(a)
a = t(scale(b))
r
pheatmap(a, cluster_rows = F, cluster_cols = F)
最后畫出來就是這樣, 明顯都是染色體兩端多
但這里有一個問題:
去重并累加
awk -F "\t" '{sum[2]+=$3}END{for(c in sum){print c,sum[c]}}' tmp3 | sort -Vk 1 > tmp4
一定要排序
累加可能不合理, 作為metaplot, 或許是每條染色體求均值, 這里先埋下伏筆
后面再改
3.前文所述問題無法解決, 所以想了想把1.中所述改一下
還是畫1.里面的圖爽室，不過不用真實的坐標和長度，全部改成染色體的相對長度
代碼就是2.所示淆攻，不過加個for循環(huán)
#Start
#獲取基因組每條染色體的長度
bioawk -c fastx '{print $name,length($seq)}' XX.Chr.fa > XX.length
python Find_A_in_B.v2.py galili.10k.fithic.tijian.density galili.length > tmp1
for i in {0..999}
do
echo $i >> tmp
done
awk '{print $1/1000}' tmp | awk '{printf "%.3f\n", int($1*1000)/1000}' > mode
rm tmp
for i in {1..26}
do
grep -w "Chr${i}" tmp1 | awk '{print $4/$2"\t"$5/$2"\t"$6}' > Chr${i}.tmp2
awk '{printf "%.3f %.3f %.3f\n", int($1*1000)/1000, int($2*1000)/1000, int($3*1000)/1000}' Chr${i}.tmp2 > Chr${i}.tmp3
awk '{print $1"\t"$1+0.001"\t"$3}' Chr${i}.tmp3 > Chr${i}.tmp4
awk -F "\t" '{sum[$1"\t"$2]+=$3}END{for(c in sum){print c,sum[c]}}' Chr${i}.tmp4 | sort -Vk 1 > Chr${i}.tmp5
python fill.gap.py Chr${i}.tmp5 mode | awk '{print $4}' | sed "1i Chr${i}" | tr "\n" "," | sed 's/,$/\n/' >> XX.tijian.everyChr.plot.csv
done
rm *tmp*
#畫圖 - 下面是R的代碼
library(pheatmap)
a = read.csv("XX.tijian.everyChr.plot.csv", header = F)
pheatmap(a, cluster_rows = F, cluster_cols = F)
畫出來如圖所示
#Find_A_in_B.v2.py
import sys
list1 = {}
with open(sys.argv[1], 'r') as f:
for line in f:
line = line.strip()
content = line.split('\t')
name = content[0]
if name in list1:
list1[name].append(line)
else:
list1[name] = [line]
with open(sys.argv[2], 'r') as f:
for line in f:
line = line.strip()
content1 = line.split('\t')
name1 = content1[0]
if name1 in list1:
for match in list1[name1]:
print(line + '\t' + match)
#fill.gap.py
import sys
list1 = {}
with open(sys.argv[1], 'r') as f:
for line in f:
line = line.strip()
content = line.split('\t')
name = content[0]
if name in list1:
list1[name].append(line)
else:
list1[name] = [line]
with open(sys.argv[2], 'r') as f:
for line in f:
line = line.strip()
content1 = line.split('\t')
name1 = content1[0]
if name1 in list1:
for match in list1[name1]:
print(line + '\t' + match)
else:
print(line + '\t' + line + '\t' + line + '\t' + "0")
