一、構(gòu)建流程:
1.UCSC參考基因組下載
2.bwa 建立index
3.samtools 生成.fai
下面以構(gòu)建人類基因組hg38為例
## 官網(wǎng)下載
cd ~/project/m6a_Mono/CHe-KK/files/genomes/hg38/bwa/
wget http://igenomes.illumina.com.s3-website-us-east-1.amazonaws.com/Homo_sapiens/UCSC/hg38/Homo_sapiens_UCSC_hg38.tar.gz
tar -xvf Homo_sapiens_UCSC_hg38.tar.gz
## 合并文件 (這里只要22條常染色體和兩條性染色體)
cd ./Homo_sapiens/UCSC/hg38/Sequence/Chromosomes
cat chr1.fa chr2.fa chr3.fa chr4.fa chr5.fa chr6.fa chr7.fa chr8.fa chr9.fa chr10.fa chr11.fa chr12.fa chr13.fa chr14.fa chr15.fa chr16.fa chr17.fa chr18.fa chr19.fa chr20.fa chr21.fa chr22.fa chrX.fa chrY.fa > hg38.fa
## 建立索引
#PS:指定算法 -a bwtsw表牢,適用于大數(shù)據(jù)窄绒;指定前綴-p hg38.fa
nohup bwa index -a bwtsw hg38.fa -p hg38.fa 2>>hg38.fa_index.log &
ls
hg38.fa hg38.fa.amb hg38.fa.ann hg38.fa.bwt hg38.fa_index.log hg38.fa.pac hg38.fa.sa #缺少hg38.fa.fai
## samtools生成 hg38.fa.fai
samtools faidx hg38.fa
上述操作構(gòu)建完成了所有的文件,可以轉(zhuǎn)移到新的文件夾中(例如~/reference/GRCh38/)初茶,便于后續(xù)使用颗祝。
二、注意事項(xiàng)
在使用的過程中恼布,可能會遇到“找不到索引文件”的報錯
[bwa_aln] fail to locate the index
[E::bwa_idx_load_from_disk] fail to locate the index files
可能有三種情況:
1.文件構(gòu)建不完整,缺少某些部分
2.電腦內(nèi)存受限制搁宾,可參考https://blog.csdn.net/weixin_40640700/article/details/116851524
3.文件名稱錯誤
例如改索引構(gòu)建后的前綴為“hg38.fa”折汞,在后續(xù)的使用中要保持一致
錯誤示例:
for f in YTHDF2.rep1 YTHDF2.rep2 YTHDF2.rep3; do
bwa aln -t 4 -q 20 ~/reference/GRCh38/hg38 ../filtering/$f.trim.c.tag.R1.fastq.gz > $f.sai
done
報錯:
[bwa_aln] 17bp reads: max_diff = 1
[bwa_aln] 20bp reads: max_diff = 2
[bwa_aln] 45bp reads: max_diff = 3
[bwa_aln] 73bp reads: max_diff = 4
[bwa_aln] 104bp reads: max_diff = 5
[bwa_aln] 137bp reads: max_diff = 6
[bwa_aln] 172bp reads: max_diff = 7
[bwa_aln] 208bp reads: max_diff = 8
[bwa_aln] 244bp reads: max_diff = 9
[bwa_aln] fail to locate the index
[bwa_sai2sam_se] fail to locate the index
正確的寫法是:
for f in YTHDF2.rep1 YTHDF2.rep2 YTHDF2.rep3; do
bwa aln -t 4 -q 20 ~/reference/GRCh38/hg38.fa ../filtering/$f.trim.c.tag.R1.fastq.gz > $f.sai
done
歡迎大家評論補(bǔ)充呀!
(每帖分享:To be whoever you wanna be!)