vcf序列提取序列及比对

使用pbsv的结果

提取序列

use strict;
use warnings;

if ($#ARGV!= 1) {
    print "Usage: perl script.pl input_file output_fasta\n";
    exit 1;
}

my $input_file = $ARGV[0];
my $output_fasta = $ARGV[1];

open(my $in_fh, '<', $input_file) or die "Can't open input file: $!";
open(my $out_fh, '>', $output_fasta) or die "Can't open output file: $!";

while (my $line = <$in_fh>) {
    chomp $line;
    my @columns = split(/\t/, $line);
    if ($#columns >= 4) {
        my $name = join("_", @columns[0..2]);
        my $sequence = $columns[4];
        print $out_fh ">$name\n$sequence\n";
    }
}

close($in_fh);
close($out_fh);

blast



#makeblastdb -parse_seqids -dbtype nucl -in xxxx.fa


#blastn -db /public/home/fengting/work/wild—38/vcf/wild-109/wild109+38/virus.taxid2lineage.all.nt.fasta -query /public/home/fengting/work/wild—38/vcf/wild-109/38ins/INS/passed/SRR13434467.fasta -out out/SRR13434467.bl -num_threads 8   -evalue 1e-5

cat ../id|while read id

do

blastn -query fasta/$id.fasta \
       -db /public/home/fengting/work/wild—38/vcf/wild-109/wild109+38/virus.taxid2lineage.all.nt.fasta \
       -out out/$id.txt \
       -outfmt 6 \
       -evalue 1e-5 \
       -num_threads 4
done

-evalue 参数
设置了期望阈值( E-value )为 1e-5, E-value 是衡量比对结果显著性的一个重要指标,它表示在随机情况下出现这样比对结果的概率,值越小说明比对结果越可靠,越不可能是随机产生的。通过设置这个阈值,可以筛选掉一些不太可靠、可能是随机比对上的结果,只保留那些更有意义的比对情况。

blast结果

blastn 是用于核酸序列比对的工具,当使用 -outfmt 6 参数时,输出结果为表格形式,每一行代表一次比对,各列的含义如下:

  1. qseqid
    含义:查询序列的标识符,即输入的待比对核酸序列的名称或编号,用于唯一标识查询序列。
    示例:如果查询序列文件中序列的标识符为 SRR13434467,则该列会显示 SRR13434467。
  2. sseqid
    含义:比对上的数据库中主题序列的标识符,代表在数据库中找到的与查询序列具有相似性的序列的名称或编号,通过该标识符可以在数据库中找到对应的完整序列信息。
    示例:若数据库中某条序列的标识符为 NC_001133.9,当查询序列与该条序列比对上时,此列会显示 NC_001133.9。
  3. pident
    含义:比对的一致性百分比,即查询序列与主题序列在比对区域内相同碱基的比例,反映了两条序列之间的相似程度,取值范围为 0 到 100。
    示例:若查询序列与主题序列在比对区域内有 80% 的碱基是相同的,则该列显示 80.00。
  4. length
    含义:比对的长度,指查询序列与主题序列比对上的区域的长度,以碱基数量表示。
    示例:如果比对区域的长度为 100 个碱基,则该列显示 100。
  5. mismatch
    含义:比对中的错配数,即在比对区域内查询序列与主题序列不相同的碱基数量。
    示例:若比对区域内有 10 个碱基不同,则该列显示 10。
  6. gapopen
    含义:比对中的空位开启数,指在比对过程中为了使两条序列更好地比对而插入的空位数量。
    示例:若在比对时插入了 2 个空位,则该列显示 2。
  7. qstart
    含义:查询序列比对上的起始位置,从查询序列的 5' 端开始计数,表明比对区域在查询序列中的起始碱基位置。
    示例:若比对从查询序列的第 10 个碱基开始,则该列显示 10。
  8. qend
    含义:查询序列比对上的终止位置,同样从查询序列的 5' 端开始计数,代表比对区域在查询序列中的结束碱基位置。
    示例:若比对到查询序列的第 100 个碱基结束,则该列显示 100。
  9. sstart
    含义:主题序列比对上的起始位置,从主题序列的 5' 端开始计数,指出比对区域在主题序列中的起始碱基位置。
    示例:若在主题序列中比对从第 20 个碱基开始,则该列显示 20。
  10. send
    含义:主题序列比对上的终止位置,从主题序列的 5' 端开始计数,确定比对区域在主题序列中的结束碱基位置。
    示例:若比对到主题序列的第 120 个碱基结束,则该列显示 120。
  11. evalue
    含义:期望值,是衡量比对结果显著性的指标,表示在随机情况下出现这样比对结果的概率。E 值越小,说明比对结果越可靠,越不可能是随机产生的。
    示例:若 E 值为 1e-10,则该列显示 1e-10,表示该比对结果随机出现的概率极低,比对结果较为可信。
  12. bitscore
    含义:比对得分,是基于比对的一致性、长度等因素计算得到的一个得分值,用于评估比对结果的质量,得分越高说明比对结果越好。
    示例:若比对得分是 100.5,则该列显示 100.5。
    通过对 blastn 输出结果中各列含义的理解,可以更好地分析和评估查询序列与数据库中序列的比对情况,从而获取有价值的生物学信息。
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,313评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,369评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 159,916评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,333评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,425评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,481评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,491评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,268评论 0 269
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,719评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,004评论 2 328
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,179评论 1 342
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,832评论 4 337
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,510评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,153评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,402评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,045评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,071评论 2 352