I am trying to use your samclip tool to keep all (PacBio HiFi) reads that map on a genome assembly with more than 1000 bp clipped (this is to try to validate a genome assembly by detecting regions where there are potential misassemblies).
Here is the way I proceeded (using samclip 0.4.0):
samclip --invert --max 1000 --ref ref.fasta < HiFimapped.sam > HiFimapped_samclip.sam
However when I check the resulting SAM I see a lot of reads that made it into the file even though they have much shorter clippings:
awk -F "\t" '{print $6}' HiFimapped_samclip.sam|head -20
1S101M1I159M1I109M1I5M1D29M1I83M1D166M1I376M1I15M1I67M1I44M1I31M1I56M1I83M1D27M1D585M1I40M1D10M1I145M2D7M1I11M1I51M1D12M1D23M1I89M2I211M1I26M1D51M1D10M1I32M1I20M1I53M1I58M1D234M1D86M1D113M1I16M1I220M1D9M1D277M1I127M1D97M1I57M1I77M1I201M1D27M1D27M1D133M1D513M1I35M1I27M1D196M2I138M1D165M1I8M1I131M2I37M1D34M1I18M1D190M1I71M1I132M1I899M1D32M1D315M1D124M1D8M1D133M1I92M1I97M1D65M1I38M1D25M1I111M1D37M1I253M1I201M1I83M1D252M1I77M1D309M1I261M1I78M1D95M1I350M1D24M1I63M1I27M1D156M1D116M1D62M1D102M1D43M1I42M1I390M1I9M1D192M1D134M1I435M1D224M1I1047M1I297M1D96M1D88M1I22M1D8M1I8M1I43M1I101M1D15M1I16M1I19M1I114M1I127M1D110M1D193M1D442M1D850M1D332M1I185M1D347M1I192M1I819M1I271M1D163M1D30M1I46M1D271M1D103M1D46M1I446M1I498M1D16M1D134M
this one seems to have a single based clipped at the beginning
1S83M2I198M1I91M1I33M1I337M1I163M1I141M1I6M1I25M1I123M1I50M1I19M1I20M1I55M1I298M1I179M1I7M1I41M1I201M1I79M1I71M2I120M1I59M1I31M1I162M1I67M1I143M1I74M1I71M1I71M1I28M1I26M2I41M1I98M3I62M2I33M1I73M1I3M1I9M1I15M1D309M1I26M1I146M1D27M1I14M2I53M1I16M2I116M1I24M1I37M2I64M1D54M1D27M1I46M1I186M1I5M1I25M1I81M1I54M1I30M1I60M3I34M1I8M1I3M1I87M1I98M1I95M2I17M1I48M1I102M1I82M1I21M1I63M1I144M1I59M1I31M1I149M1I236M1I8M1D81M1I67M1I75M1I75M1I4M1I42M2I89M1I150M1I16M1I138M1I13M1I9M2I54M2I302M1I26M1I18M1I15M1I34M2I46M1I122M1I68M1I103M1I307M1I232M1I29M1I27M1I130M1I93M1I158M1I64M1I63M1I125M1I32M1I87M1I210M2I273M2I68M1I140M1I101M3I936M1I8M1I53M1I200M1D87M1I13M1I26M1I20M2I93M1I178M1I226M1I14M1I94M1I40M1I171M1I273M1I44M1I814M1I12M1D10M1I28M1I245M1I41M1I38M1I74M1D66M2I25M1I77M1I180M2I40M1I125M1D196M1I42M1I406M1I5M1I60M2I36M1I9M1I31M1I20M1I21M1I9M1I71M1I74M1I56M1I158M1I14M1I27M1I36M
this one too
179M1I1232M1I373M1D66M1I167M1D351M1I52M1I117M1I260M1I377M1I138M1D315M1D17M1I50M1D732M1D158M1I159M1I66M1I153M1I63M1D9M1I1768M1I90M1D398M1I515M1I110M1I45M1I754M1D276M1I137M1D494M1I285M1I60M1I148M1I24M1I734M1D147M1I377M1I352M1I700M1D16M1D1172M1I219M1D66M1I386M1I182M1I491M1I83M1I856M1D167M1I30M1D321M1I335M1I633M1D16M1I57M1I147M1D171M1I102M1I37M1D427M1I129M1D149M1D35M1I187M1I17M1I163M1D187M1I66M1I190M1I78M1I233M1I215M1I121M1I63M7S
this one has only 7 bases clipped at the end
Hi Torsten,
I am trying to use your samclip tool to keep all (PacBio HiFi) reads that map on a genome assembly with more than 1000 bp clipped (this is to try to validate a genome assembly by detecting regions where there are potential misassemblies).
Here is the way I proceeded (using samclip 0.4.0): samclip --invert --max 1000 --ref ref.fasta < HiFimapped.sam > HiFimapped_samclip.sam However when I check the resulting SAM I see a lot of reads that made it into the file even though they have much shorter clippings: awk -F "\t" '{print $6}' HiFimapped_samclip.sam|head -20
1S101M1I159M1I109M1I5M1D29M1I83M1D166M1I376M1I15M1I67M1I44M1I31M1I56M1I83M1D27M1D585M1I40M1D10M1I145M2D7M1I11M1I51M1D12M1D23M1I89M2I211M1I26M1D51M1D10M1I32M1I20M1I53M1I58M1D234M1D86M1D113M1I16M1I220M1D9M1D277M1I127M1D97M1I57M1I77M1I201M1D27M1D27M1D133M1D513M1I35M1I27M1D196M2I138M1D165M1I8M1I131M2I37M1D34M1I18M1D190M1I71M1I132M1I899M1D32M1D315M1D124M1D8M1D133M1I92M1I97M1D65M1I38M1D25M1I111M1D37M1I253M1I201M1I83M1D252M1I77M1D309M1I261M1I78M1D95M1I350M1D24M1I63M1I27M1D156M1D116M1D62M1D102M1D43M1I42M1I390M1I9M1D192M1D134M1I435M1D224M1I1047M1I297M1D96M1D88M1I22M1D8M1I8M1I43M1I101M1D15M1I16M1I19M1I114M1I127M1D110M1D193M1D442M1D850M1D332M1I185M1D347M1I192M1I819M1I271M1D163M1D30M1I46M1D271M1D103M1D46M1I446M1I498M1D16M1D134M this one seems to have a single based clipped at the beginning
1S83M2I198M1I91M1I33M1I337M1I163M1I141M1I6M1I25M1I123M1I50M1I19M1I20M1I55M1I298M1I179M1I7M1I41M1I201M1I79M1I71M2I120M1I59M1I31M1I162M1I67M1I143M1I74M1I71M1I71M1I28M1I26M2I41M1I98M3I62M2I33M1I73M1I3M1I9M1I15M1D309M1I26M1I146M1D27M1I14M2I53M1I16M2I116M1I24M1I37M2I64M1D54M1D27M1I46M1I186M1I5M1I25M1I81M1I54M1I30M1I60M3I34M1I8M1I3M1I87M1I98M1I95M2I17M1I48M1I102M1I82M1I21M1I63M1I144M1I59M1I31M1I149M1I236M1I8M1D81M1I67M1I75M1I75M1I4M1I42M2I89M1I150M1I16M1I138M1I13M1I9M2I54M2I302M1I26M1I18M1I15M1I34M2I46M1I122M1I68M1I103M1I307M1I232M1I29M1I27M1I130M1I93M1I158M1I64M1I63M1I125M1I32M1I87M1I210M2I273M2I68M1I140M1I101M3I936M1I8M1I53M1I200M1D87M1I13M1I26M1I20M2I93M1I178M1I226M1I14M1I94M1I40M1I171M1I273M1I44M1I814M1I12M1D10M1I28M1I245M1I41M1I38M1I74M1D66M2I25M1I77M1I180M2I40M1I125M1D196M1I42M1I406M1I5M1I60M2I36M1I9M1I31M1I20M1I21M1I9M1I71M1I74M1I56M1I158M1I14M1I27M1I36M this one too
179M1I1232M1I373M1D66M1I167M1D351M1I52M1I117M1I260M1I377M1I138M1D315M1D17M1I50M1D732M1D158M1I159M1I66M1I153M1I63M1D9M1I1768M1I90M1D398M1I515M1I110M1I45M1I754M1D276M1I137M1D494M1I285M1I60M1I148M1I24M1I734M1D147M1I377M1I352M1I700M1D16M1D1172M1I219M1D66M1I386M1I182M1I491M1I83M1I856M1D167M1I30M1D321M1I335M1I633M1D16M1I57M1I147M1D171M1I102M1I37M1D427M1I129M1D149M1D35M1I187M1I17M1I163M1D187M1I66M1I190M1I78M1I233M1I215M1I121M1I63M7S this one has only 7 bases clipped at the end
Am I doing something wrong?
Thanks a lot in advance for your help!