Question: Hello folk, I'm beginner in perl and I'm trying to extract the whole column with a set of specific patterns in tab-delimited file (very large). I saw some scripts but only when you know which column are you interested in (small file).

My file structure looks like this:

##info1
##info2
##info3
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	ID01	ID02	ID03	etc...
3	66894	rs9681213	0	1	.	PASS	.	GT	0|1	0|1	0|1	etc...
3	95973	rs1400176	0	1	.	PASS	.	GT	1|1	1|1	1|1	etc...
3	104972	rs990284	0	1	.	PASS	.	GT	0|1	0|1	0|0	etc...
3	114133	rs954824	0	1	.	PASS	.	GT	1|1	1|1	1|1	etc...
and so on...

Let's suppose that I want to extract four columns (ID01, ID03, ID45, ID80). These IDs could be (optionally) in another text file for parsing. I suppose I have to storage these patterns in an array, right!?

I tried to split the @_ data (for getting the lines) but I couldn't get the correct column $_[?] to print off. The file is too big and I can't see in which column number they are. I really get stuck...

That's the code I've been working on:

#! usr/bin/perl
use strict;
use warnings;

my $file = $ARGV[0] || die;
my @fname = split(/_/, $file);	#my file name is chr3_etc.vcf

open(FILE, '<', "$file") || die;
open(OUT,'>', "$fname[0]_NAD.vcf");

my @NAD_ids = qw(NAD15 NAD54 NAD55 NAD56 NAD57 NAD58 NAD59 NAD64 NAD93 NAD93 NAD98);

while () {
next if m/##/;	#exclude the additional information above the table
chomp; @_=split;	#getting the columns
if (@NAD_ids) {
print OUT "$_[0]\t$_[1]\t$_[2]\t$_[3]\t$_[4]\t$_[5]\t$_[6]\t$_[7]\t$_[8]\t$_[($NAD_ids[0])]" . "\n";	# see the last one (problem);
}

}

close FILE; close OUT;
exit;

Please consider the "NADs" like the IDs.

Answer: You are very close and I appreciate the effort. I'll share the thing you are missing and some other helpful tips because you basically have all the logic right. (I personally think it would be better to use vcf-tools or one of the associated VCF modules to work on the indexed files more efficiently. That being said, let's work with what you posted.) The most important thing to remember is to use lexical variables to limit the scope of your actions. When you are using loops and modifying the topic `$_` or trying to directly modify things like `@_` (this special notation only applies to subroutine arguments, so it was misplaced) it can get confusing. That is one reason it's best to work with lexical variables and avoid the magic (also, readability). Exceptions are in map/grep/sort blocks. Also important is that bare filehandles have global scope, so they are to be avoided in all cases. These can cause problems from a large distance in the code.

For the question, it sounds like you want to pick genotypes from one file and select those from a VCF file, right? It's not completely obvious at first because your examples don't match up. Also, there were duplicate IDs in your example array, which may or may not matter but it's important to check your inputs (you can solve this problem using the same module I used below using the `uniq` method). I modified the files to demonstrate, otherwise there was nothing to be done. The ID file:

https://gist.github.com/sestaton/7426cb038c5ae09c7ae6853735d597d8

And the VCF:

https://gist.github.com/sestaton/5d527b0b5accdba812d1267d971a9a40

Here is some code:

https://gist.github.com/sestaton/d9de56fe1b3b0e78208ff515cbe796d8

The first thing to do is read in the IDs, which is what the `map` block is doing. What that is doing specifically is removing newlines and spaces so the IDs will match those in the file, and then assigning the keys a value of 1 (to say it's 'defined' in Perl-speak). The space/newline issue is important in any programming language. The rest is basically what you had with some minor changes. To pick the correct genotypes you find the column (array index), then test if any matches were found with `if (@indexes) {`, and if not you stop reading to be efficient (literally, you say that's the `last` line we need). The data is printed with an "array slice" and `join`, and these two idioms will save you a lot of typing. If you are unsure about what is happening here, use [Data::Printer](https://metacpan.org/pod/Data::Printer) (the usage of which will be very familiar if you know Ruby) or [Data::Dump::Color](https://metacpan.org/pod/Data::Dump::Color) (personal preference) to peek and the data structures in the script.

Running this prints:

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NAD93	NAD98	NAD58
3	66894	rs9681213	0	1	.	PASS	.	GT	0|1	0|1	0|1
3	95973	rs1400176	0	1	.	PASS	.	GT	1|1	1|1	1|1
3	104972	rs990284	0	1	.	PASS	.	GT	0|1	0|0	0|1
3	114133	rs954824	0	1	.	PASS	.	GT	1|1	1|1	1|1

That's what I think you wanted. If that's not exactly what you were asking then I'm confident you can modify it to your needs. That will likely be good enough, but as I mentioned above, I'd personally use a toolkit or some kind of VCF library for this since it has to be a solved problem (I recall doing this recently and I’ll update if I remember the command).