How to save the results of norm to txt
1
0
Entering edit mode
8 months ago
Serij´s • 0

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 • 1.2k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode
8 months ago
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
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

What $SHELL are you using?

ADD REPLY
0
Entering edit mode

I am using Ubuntu but with zsh shell

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

yes, unless OP wants to save the vcf

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 1134 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6