0
$\begingroup$

I'm trying to study fetal fraction by using a seqff script from here. The instruction says that I need a headless sam file to work so after alignment with bwa, I process my sam file by using samtools:

samtools view -S -q 20 -F 4 target.sam > treated.sam

This process also remove the @SG header since there is no reference file. I also use grep to remove header of the target.sam got after alignment for comparison. However, when I compare, I noticed that ENet calculation on the treated.sam returned with NA results, so half of the method to calculate fetal fraction of this script was gone. I have tested it in multiple cases and it seems that without using samtools, the script worked well but with samtools Enet calculaton not worked. Can someone tell me what is wrong with the script? I'm using R 4.0 and all modification I made to the script is on the args since it cannot recognize args I input with anymore.

p/s: this is few lines of headless target.sam made by using grep

W2003120038L1C001R0050000009    0   chr10   6490952 0   50M *   0   0   GAGCAGAAAAACTGGAAACTCTAAAAAGCAGAGTGCCTCTCCTCCTCCAA  IBGIBAB<ABBCA9IB@BHDIA=ABBB8IB9A;DGII@I8IIAII<IDCA  XT:A:R  NM:i:0  X0:i:151    XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000037    0   chr5    172722817   37  50M *   0   0   GTCCTAAGACAAAGACCTCAGTGCTTATAATGTAAGTGGTCACTGGGCTG  ICIICBB:BIBBBIAICDIAIAIICDBDAABIDBBICEIAIAFDH?IIAI  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7T42
W2003120038L1C001R0050000062    0   chr5    84724587    37  50M *   0   0   TAAGAATTGTGCATAACCAGTAGCCCATGTTATCTTGCAGAAGCATTTGT  ABBI@BDDIDIIADABHEAIC@IIIBBAID@ADIDCA:BC=BIIAAAAIC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000067    16  chr5    45648067    37  50M *   0   0   TTTATTTAAATACAAGCCAATTAAATGATTTTATTTTACATTTTCTAACC  AABAAABCAAABH@CFHIDDB@BCDAID@BBACBB@BBIDB@2BIBC@II  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000184    16  chr4    103629659   37  50M *   0   0   CAGTGCTAAGAGGGAAGTTTATAGCAACAGATACTTACATCCAAAAAGTA  FDFBIIADAFDICICCIAABDBAFIADIDFDABIBA@I@CI3?DCCDIBC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000211    0   chr3    157734199   37  50M *   0   0   TCTGAAATTGATCTTTATCCCACATCCACGTATCCTAAAATAATAAAACA  AFDBBBB*CIBDFACABDIIBB<1CI=BFIDBAFEABBAADBBABAABIB  XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:7G15T26

And this is from treated.sam, the sam file be processed by samtools

W2003120038L1C001R0050000037    0   chr5    172722817   37  50M *   0   0   GTCCTAAGACAAAGACCTCAGTGCTTATAATGTAAGTGGTCACTGGGCTG  ICIICBB:BIBBBIAICDIAIAIICDBDAABIDBBICEIAIAFDH?IIAI  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7T42
W2003120038L1C001R0050000062    0   chr5    84724587    37  50M *   0   0   TAAGAATTGTGCATAACCAGTAGCCCATGTTATCTTGCAGAAGCATTTGT  ABBI@BDDIDIIADABHEAIC@IIIBBAID@ADIDCA:BC=BIIAAAAIC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000067    16  chr5    45648067    37  50M *   0   0   TTTATTTAAATACAAGCCAATTAAATGATTTTATTTTACATTTTCTAACC  AABAAABCAAABH@CFHIDDB@BCDAID@BBACBB@BBIDB@2BIBC@II  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000184    16  chr4    103629659   37  50M *   0   0   CAGTGCTAAGAGGGAAGTTTATAGCAACAGATACTTACATCCAAAAAGTA  FDFBIIADAFDICICCIAABDBAFIADIDFDABIBA@I@CI3?DCCDIBC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000211    0   chr3    157734199   37  50M *   0   0   TCTGAAATTGATCTTTATCCCACATCCACGTATCCTAAAATAATAAAACA  AFDBBBB*CIBDFACABDIIBB<1CI=BFIDBAFEABBAADBBABAABIB  XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:7G15T26
W2003120038L1C001R0050000212    0   chr4    157022243   37  50M *   0   0   AAAAATAAAAAAGCATTCTGCCCAGATCTCTCAACATTATGTTATTTCTT  4<AA?D2,AAA@IF@DDIDIIDIA2AD=DBAF@BIADA<DICD@DADFD@  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000258    0   chr13   106093727   37  50M *   0   0   AATGGTTACTTTATGTTTATATTTCAACAGGGAATATTTGAAATATGGCA  97?IIDA*FDAD<@H@AD<@<?AAG8??,I;I@:DADDAG7@@7(;IH6?  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7C42
$\endgroup$
5
  • $\begingroup$ @gringer I have added few lines of both file on my post. Also I have to perform quality filter to filter out bad alignments results $\endgroup$
    – khanhlpbao
    Apr 22 at 7:35
  • $\begingroup$ If the grep one works, and the non-grep one doesn't, my first step would be to exclude any confounding factors (like quality filtering) to try to get them both to look the same (which they don't, at the moment). Is there any improvement if you don't include the quality filter (either in the files looking similar, or the script working better)? $\endgroup$
    – gringer
    Apr 22 at 11:26
  • $\begingroup$ If you think there is a formatting problem, you might need to put small subsets of the samtools file through until you can identify the problematic line. $\endgroup$
    – swbarnes2
    Apr 22 at 23:11
  • $\begingroup$ @swbarnes2 thanks, do you know any command or program that can check for the different? Because the each sam file is more than 4GB, finding what is different by compare with text editor is impossible to say at least $\endgroup$
    – khanhlpbao
    Apr 23 at 4:01
  • $\begingroup$ One you figure out the problem in one file, you should be able to fix it for the others. $\endgroup$
    – swbarnes2
    Apr 23 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.