根据SNP位置提取附近序列，再blast，进行不同基因组位置的替换
两个位置之间大点，太小blast不出来，bed文件的染色体标注在不同参考基因组不一样
vi 2.bed
Chr1 7720000 7720200
Chr1 275000000 275000200
Chr1 92800000 92800200
Chr3 162000000 162000200
Chr1 109800000 109800200
Chr2 13900000 13900200
Chr7 14500000 14500200
bedtools
bedtools getfasta -fi /data4/husk_RNA-seq_wsk/HAPMAP3/B73_RefGen_v2.fa -bed 2.bed -fo hn_v2.fa
Chr1:7720000-7720200
GGTGACTAATTAAGCATGCGGTCAATGTAGCTTGGACACTAGTAAACTAGACAAAAAATAATAATATTTGTGACACAAGAATCAGTTGAATATTGCAAATTCACATAAAGCTGGTCCAATGTGGACATAAGCTA>
Chr1:275000000-275000200
aagagtatttttccctctagttttaatggaaaaattaattccgaggaggatcaaagatgttatcgcttaaactatctagctcttgacatcttaactaattctcttagcaaagaagattatcgtgccttcatatc>
Chr1:92800000-92800200
ctgaaatttctggtaccgaagctccagagcaaaaaactgaagccgaagctgggccttcagagcccgcaaagataaaatccttggagagtgaggtggaaaaaataacaaagccaacttttgttgaagaaatcggt>
Chr3:162000000-162000200
atgtcatccttccccctgggggagatcatccagtgccgaggggcctcgggcaggatcgcaaagtgggcagtggaaatcatgggcgaaacgatctcgttcgcccctcggaaggccatcaagtcccaggtgttggc>
Chr1:109800000-109800200
aagtgctttctagtgaattatcttcttgtcatgaaaacattgctactttaaagagtgttaatgatgacttaaatgctaaactagaagtagctagtaaatcaacatcttgtgtagaaattgttgaaacgtgcact>
Chr2:13900000-13900200
tctagtagttaaagaatcaagatctattcctatccgtggaagaacctcaagttttcaataagcaaggcaagtggacttctccctctgcatactctgtttaatcccaaacaatacatttaataaagtttattctt>
对参考基因组建库
makeblastdb -dbtype nucl -parse_seqids -in ../Zea_mays.B73_RefGen_v4.dna.toplevel.fa -out B73_database
本地blast,很快
blastn -query hn_v2.fa -out seq.blast -db B73_database -outfmt 6 -max_hsps 1 -num_threads 8
保留每条序列，最高的一条，检查两个版本位置是否一样？
awk '!a[$1]++{print}' seq.blast > uniq.blast
B73_V1 参考基因组有问题，需要整理一下再构建索引
error1
seqkit subseq --bed 5.bed -o 5.fasta /data4/husk_RNA-seq_wsk/ZmB73_AGPv1.fa
[INFO] read BED file ...
[INFO] 24 BED features loaded
[INFO] create or read FASTA index ...
[INFO] read FASTA index from /data4/husk_RNA-seq_wsk/ZmB73_AGPv1.fa.seqkit.fai
[WARN] 0 records loaded from /data4/husk_RNA-seq_wsk/ZmB73_AGPv1.fa.seqkit.fai, please check if it matches the fasta file, or switch on the flag -U/--update-faidx
error2
different line length in sequence: chr1. Please format the file with 'seqkit seq'
将此文件fasta序列转换成100个碱基一行输出
seqkit seq -w 100 ZmB73_AGPv1.fa > ZmB73_AGPv1_2seq.fa
seqkit subseq --bed 5.bed -o 5.fasta /data4/husk_RNA-seq_wsk/ZmB73_AGPv1_2seq.fa
