Entering edit mode
8.2 years ago
beausoleilmo
▴
600
I want to use samtools to create a .bcf file using this command:
samtools mpileup -C 50 -E -t SP -t DP -u -I -f /my/path/Geospiza_fortis.scaf.noBacterial.fa -b bam_list.txt > output_test.bcf
But I get this error:
[fai_fetch_seq] The sequence "JH739887" not found
I know that the problem is that the name "JH739887" doesn't exist in the reference genome. But How I'm I supposed to change this name so that it can be mapped on the reference genome?
head Geospiza_fortis.scaf.noBacterial.fa
>Scaffold410 275
TGCATTAATATGAGTGTGTGCTGCAAAAGTTCAGGTCATGGTCCGATCATACTTCACATTTTGGTAGCACTTTAAGCAGAGATCGGTTATCCCATTCTGTGGAAGACTCAACACTATCATAAGGTCCCACAGTTTTATTATCCCTCTGCCTCCCGGAATGCCCCCGGCAGTGAGGGGTACCATCTTCTCAGCAGTAAGGATATTCTTCAGGAGTTCCGTGTGAGCTTTCCCGGATTTAGTTCCATTTTTTAAATACTTCCCAATTCTTTGCTTTG
>Scaffold430 374
CTTTGTTAACTGAAAGAGCCTCTAAGTAGATGACCAGTGCTCAGTTAGTACAGTATGAATTTTGTTTAATGGAACAGGAAGATTTAGTATTGAGAAGCGGTTAAGGGTTTAACCCAGCCTCCTGTCTGAATGGACCTGAAGAGGGGGGCCGGGAAGAAACCCATGACTGCATTAAAGTGATAGATCTCCAGACATGGGCTAGGGAAGATTTACAAGACACTCCCTGGCCTGAGGGAGAAAATATGTTTATTGATGAGTCTTCAAGGGTGGCAGAAGGGAAGCGATTTACAGGATACACAATCATTAATGGAAGGAAATTAAAGGAAGGGGGGAGATTGTCACCCACCTGGTCAGTTCAGACAGCAGAGCTGTAT
head Geospiza_fortis.scaf.noBacterial.fa.fai
# Scaffold410 275 17 275 276
# Scaffold430 374 310 374 375
# Scaffold1010 597 703 597 598
# scaffold260 820 1314 820 821
# scaffold980 554 2148 554 555
# scaffold440 1463 2716 1463 1464
samtools view Bamboum_sorted.bam | head
24016062 16 JH739887 2197 37 94M * 0 0 ACCATTTAAAACACTCTTCAGGGAAGGCTTTAGCTTTTCATTCAGTTTTGTATTTTACTACATGAATATAAAGCAGTTTTGCTTCTGTGTTTTT >@BB@>>>;;?>;@A>AAABBBBBBBBBBBBBBBBBBA?A?9<===<<BAABBABA=BBABBBBBBBBCBB>CCCACCC73?@2+<2+0+A<+= NH:i:1 NM:i:3 RG:Z:JP3750_for_sorted
27522854 0 JH739887 2921 37 94M * 0 0 AATTTCCAACACTCTGCAAGAAGCTGCAAGAGTTCGCCTGCACCTGTGTCTTATTGCCACAAGAGCTCTCTGTTCTATACTCTGAAATTCTCAC BCFFFFFHHHHHJJJJJJJJIJIJIJJJJJJJFGHIIJJJIJJJIIIIIHHIJJJJJJJJJJHIJJJJHHHHHHFFFFFFFEEEEEEEDDEDED NH:i:1 NM:i:0 RG:Z:JP3923_for_sorted
31888152 0 JH739887 4964 37 94M * 0 0 AAAGCCCACCTGCATTATACAACATGAAGCCGTAGAATAGACAGTCTGATGGTTTACAATGTTTAGGAGATGTTACCATGTTAAGGATCATTAT C@FFDFFHHHHHJIJJJJJJJJJJJIJJJJJJGIJJJJJ*?GHIEHIJIHIIHIIJJIJJJHII===CHEHGHHFFFFFFEEEDEEEDDDDEEE NH:i:1 NM:i:5 RG:Z:JP3752_for_sorted
Had you aligned against the fasta file that you're giving to samtools you wouldn't have these issues.
@beausoleilmo: Heed the advice of someone who knows best :)
Someone sent me this file containing the different possible names. You can use a script to change the name of the files in either of the files (FASTA __or__ BAM)
ask your collaborator the REFerence fasta file he used with BWA. period.