在GEO數(shù)據(jù)庫下載SHARE-seq的數(shù)據(jù)后菠净,發(fā)現(xiàn)這些fragment文件與CellRanger輸出的fragment文件存在一下幾點不同:
- SHARE-seq的fragment文件只有4列,并且存在重復行屈梁,重復的行數(shù)為該條read的counts數(shù)
- 沒有以.tbi結尾的索引文件
- 所提供的文件壓縮格式不是BGZF嗤练,這樣會導致在建立索引時報錯
- 行是亂序的,后面行的起始位點小于前面行的起始位點在讶,這樣會導致在建立索引時報錯
- 細胞名稱與所提供矩陣中的細胞名稱不同煞抬,fragments文件中為
,
連接,而矩陣中為.
連接
以下為處理流程:
- 首先安裝htslib构哺,需要使用其中的tabix工具對fragment文件進行.tbi索引文件創(chuàng)建
# 下載革答,手動編譯安裝
wget https://github.com/samtools/htslib/releases/download/1.14/htslib-1.14.tar.bz2
tar -jxvf htslib-1.14.tar.bz2
cd htslib-1.14
./configure --prefix=/local/txm/software/htslib
make
sudo make install
# 添加環(huán)境變量
vi ~/.bashrc
export PATH=/local/txm/software/htslib/bin:$PATH
source ~/.bashrc
- 針對以上幾點問題,首先對所提供的fragment.bed.gz文件解壓曙强,將
,
全部替換為.
gunzip GSM4156599_brain.atac.fragments.bed.gz
sed -i "s/,/./g" GSM4156599_brain.atac.fragments.bed
然后對GSM4156599_brain.atac.fragments.bed文件依次進行以下操作:
1 . 對文件按照前三列進行排序
2. 去除重復行残拐,并統(tǒng)計重復次數(shù)(該read的counts數(shù))
3. 默認統(tǒng)計的次數(shù)在第一列,將其放至最后一列(標準fragment文件格式)碟嘴,并使用制表符隔開
4.最后使用bgzip進行壓縮命令如下
sort -k1,1V -k2,2n -k3,3n GSM4156599_brain.atac.fragments.bed | uniq -c | awk '{print $2,$3,$4,$5,$1}' OFS="\t" | bgzip -c > GSM4156599_brain.atac.fragments.mybed.gz
- 最后對mybed.gz文件進行索引文件生成
tabix -0 -p bed GSM4156599_brain.atac.fragments.mybed.gz
大功告成溪食,在下游分析中記得再把rna中的細胞名全部換成atac的就可以保證細胞名是統(tǒng)一的了
20230201 Bing Ren scATAC-seq human atlas fragments文件處理記錄
### ADRUQ
# 首先解壓原來的fragments文件,因為不是BGZF格式
gunzip GSM5589391_pancreas_SM-ADRUQ_rep1_fragments.bed.gz
# 依次排序娜扇、去掉最后一列(第六列)错沃、壓縮
sort -k1,1V -k2,2n -k3,3n GSM5589391_pancreas_SM-ADRUQ_rep1_fragments.bed | awk '{print $1,$2,$3,$4,$5}' OFS="\t" | bgzip -c > GSM5589391_pancreas_SM-ADRUQ_rep1_fragments.mybed.gz
# tabix生成索引
tabix -0 -p bed GSM5589391_pancreas_SM-ADRUQ_rep1_fragments.mybed.gz
參考
http://www.htslib.org/doc/tabix.html
http://www.reibang.com/p/912c81d71045
https://blog.csdn.net/biocity/article/details/83274985