This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to save the results of norm to txt

Hello,

I already asked the question the comments of my last question : link to last question

But I think the question warrents its own post.

I am using the following command to detect missing values in REF:

bcftools norm --check-ref w --fasta-ref in.fa

I want to save the result of the analysis in a txt-file. I want to creat a loop later which takes multiple gvcf files and detects missing values in every REF-Column. The results of all files should be than exported to one txt-file.

I grepl the last line containing the analysis-results using grep:

grep "Lines total/split/realigned/skipped"

and add the output to a txt-file instead of printing it in terminal:

bcftools norm --check-ref w --fasta-ref in.fa input.g.vcf.gz | grep "Lines total/split/realigned/skipped" >> missing_data.txt

I also tried tree, > and 2>&1.

Nothing worked.

Is there a better way to 1. grep the last line and 2. export the analysis-results to txt.

bcftools variant-calling vcf
• 2,716 views • link
updated by Ram Ram • written by Serij´s • Serij´s

If you already asked it, why ask a new question instead of just waiting for a response?

• 0 views • link
written by Ram Ram

1 answer 

bcftools norm --check-ref w --fasta-ref in.fa-O z  -o output.vcf.gz  input.g.vcf.gz  2>&1 | grep "Lines total/split/realigned/skipped" >> missing_data.txt
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thank you for the quick help.

Sadly, it also doesnt seem to work for me. Maybe I did it wrong or I am missing something, but right now the missing_data.txt is being generated, the code runs through but the missing_data.txt remains empty after execution

bcftools norm --check-ref w --fasta-ref "/media/x/My Book/x_data/data/Homo_sapiens_assembly38.fasta" -O z -o "output.vcf.gz"  "input.g.vcf.gz" 2>&1 | grep "Lines total/split/realigned/skipped" >> "/home/X/Desktop/missing_data.txt"

I was also wondering, if there was a way to execute the code without generating new vcf files, since my input data is really large.

• 0 views • link
written by Serij´s • Serij´s

Why not just set -o to /dev/null?

• 0 views • link
written by Ram Ram

I tried that. That didnt work for me.

using script command worked, but its not an elegant solution.

output_line=$(script "$temp_file" -c "bcftools norm -f $fasta_ref $gvcf_file | grep 'Lines total/split/realigned/skipped'")

# Append the output line to missing_data.txt
echo "$gvcf_file" >> "$output_file"
echo "$output_line"| sed -n '2p' >> "$output_file"

If there is a better way to export only the following part of the resulst, that would be great

Lines   total/split/realigned/skipped:  453061507/0/0/0
• 0 views • link
written by Serij´s • Serij´s

works on my machine:

bcftools norm -f rotavirus_rf.fa --check-ref w rotavirus_rf.vcf.gz -o /dev/null  2>&1 | grep 'total/split/realigned/skipped:' > out.txt && cat out.txt
Lines   total/split/realigned/skipped:  45/0/5/0
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

What $SHELL are you using?

• 0 views • link
written by Ram Ram

I am using Ubuntu but with zsh shell

• 0 views • link
written by Serij´s • Serij´s

zsh can do odd things with redirects and I'm not savvy enough to debug that remotely, so try switching to bash and then running the command given by Pierre.

• 0 views • link
written by Ram Ram

I didnt know that, since >> and tree usually worked under zsh. But thank you for the tip. I will try it out and let you know if it worked or not.

• 0 views • link
written by Serij´s • Serij´s

ok, Thank you. I tried it and it worked in bash but not zsh.

• 0 views • link
written by Serij´s • Serij´s

There is a setting that interferes with 2>&1 when you use a pipe on zsh. This is from my vague recollection, so I might be off on specifics.

EDIT: Found the keyword, it's called MULTIOS. Read up on it when you can: https://unix.stackexchange.com/a/672946/135331

• 0 views • link
written by Ram Ram

No. I have 700 gvcf files with half a billion rows in each gvcf. Thats why I only want to save the results of the analysis and not the content of each gvcf-file into one txt.

• 0 views • link
written by Serij´s • Serij´s

Log in to answer this question.

More posts like this