二代数据序列QC,BWA,IGV可视化快速流程

1. Reads QC

使用fastqc 软件

参数：

• -t 线程数
• -o 输出文件夹名
• 最后输入要质控的序列名字

fastqc -o qc -t 15 S11.fastq# 

使用prinseq-lit软件

去除左端10bp，右端5bp，最低长度1000，输出到text.qc.fq中

参数:

• -verbose 动态输入过程
• -fastq 输入文件格式和名字
• -trim_left 左端删除多少bp
• -trim_right 右端输出多少bp
• -min_len 最短的长度
• -out_good 输出符合要求的good序列名字

perl /home/dir/mixassembly_congig/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fastq test2.fq -trim_left 10 -trim_right 5 -min_len 1000 -out_good text.qc.fq

不输出bad seq

perl /homes/dir/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fastq S${i}.fastq -trim_left 10 -trim_right 1 -min_len 1000 -out_good S${i}_qc.fq -out_bad null 

• /dir 为prinseq-lite目录

使用awk函数过滤长度

长度过滤还有可以用awk函数，一行代码即可搞定

1. 单行fasta：
   假设过滤40bp的序列

awk '{if($0~/>/) geneName=$0;if(length($0)>40) print geneName"\n"$0;}' seq

2. 多行fasta：
   假设过滤70bp的序列

awk 'BEGIN{OFS=FS="\t"}{if($0~/>/) name=$0 ;else seq[name]=seq[name]$0;}END{for(i in seq) {if(length(seq[i])>70) print i"\n"seq[i]}}' test.fasta

2. 使用BWA比对

bwa mem 比对

首先对参考基因组做索引

bwa index reference.fa

正式比对

参数

• -t 线程
  
• 后面跟read1和read2

bwa mem -t 15 refernce.fa reads_1.fq reads_2.fq >Sample.sam

使用samtools转换sam为bam

1. 使用sam view 把sam转换成bam

samtools view -bS sample.sam >sample.bam

2. 建立比对文件排序：samtools sort，生成sample.sort.bam文件
   samtools sort 输入的bam文件 sort后的bam文件名称(会自动加.bam后缀)
   
   参数：

   -@ 线程数
   -m 使用最大内存数

samtools sort sample.bam sample.sort #produce sample.sort.bam
samtools sort -@ 15 -m 15G sample.bam sample.sort

3. 对sort后的bam文件索引

samtools index sample.sort.bam #produce sample.sort.bam.bai

3. 使用IGV可视化mapping的结果

首先在genome标签里导入参考序列
然后在file里面导入sort后的bam，同时bam文件夹里面必须有bam文件的索引bai文件
1.genome import reference.fa
2.file import sample.sort.bam #must have bai file in the same folder

4. SNP Calling

使用samtools

1. 首先建立参考序列的索引，使用samtools faidx命令

samtools faidx KPN_target_genes_online_UP.fa

2. 使用mpileup生成每个位点碱基累计数据，并用view转换bcf为vcf文件
   
   参数

-g Compute genotype likelihoods and output them in the binary call format (BCF)(计算genotype likelihoods并以bcf格式输出)

-f The faidx-indexed reference file in the FASTA format(samtools faidx生成的索引.fai文件，用其他软件建索引也可以的)

-u Generate uncompressed VCF/BCF output(输出不压缩的bcf文件)

samtools mpileup -g -f KPN_target_genes_online_UP.fa s1.sosrt.bam.bam >s1.raw.bcf
bcftools view -cvg s2.raw.bcf >s2_var.vcf

使用freeways

/homes/xiaohuizou/tools/freebays/freebayes/bin/freebayes -f ref.fa s1.sort.bam >freebays/s1_freebays.vcf

用scapel 检测indel突变

for i in {13,17,18,26,32,34,37,38,39,40,47,48,61,62,66,67,70,71,72,76,83,86,92,96};do scalpeld --single --bam ../s${i}.sort.bam --bed ../jsref.bed --ref ../jsref.fa --verbose --numprocs 15 --dir s${i};done

5. SNP过滤

1. 删除特定序列的结果：

删除含有kpc的行*

g/kpc2/d # vim删除kpc2的行：
sed -i '/kpc/d' s[1-9]*.var.vcf # sed匹配模式  
sed -i '/KPC/d' s[1-9]*.var.vcf #sed直接删除

2. 过滤深度

过滤depth超过不低于30的snp

awk -F '[\t;]' '{split($8,DP,"=");if(DP[2]>=30) print $0}' s1.var.vcf > s1.var.dep30.vcf

3. 过滤info DP4信息：

• 匹配含有DP4的lie

awk -F '[\t;]' '{for(i=1;i<=NF;i++)if($i~/DP4/) print $i}' s1.var.dep30.vcf
awk -F '[\t;]' '{for(i=1;i<=NF;i++)if($i~/DP4/) {print $i;}}' s1.var.dep30.vcf

结果：

DP4=270,201,285,3

DP4=284,221,302,9

DP4=198,37,260,270

DP4=0,0,108,201

DP4=100,342,6,247

DP4=117,313,5,206

DP4=111,152,14,313

• 使用split拆分字DP4字符串

awk -F '[\t;]' '{for(i=1;i<=NF;i++)if($i~/DP4/) {split($i,ac,",")}print ac[1]}' s1.var.dep30.vcf 

结果

DP4=270

DP4=284

DP4=198

DP4=0

• 打印正负链均大于5个reads支持的位点

awk -F '[\t;]' '{for(i=1;i<=NF;i++)if($i~/DP4/) {split($i,ac,",")}if(ac[3]>=5||ac[4]>=5)print $0}' s1.var.dep30.vcf >s1.fin.vcf

加上vcf 的header：

grep '#' s1.var.vcf >title<br>
awk '{if(FNR==1)print system("cat title")$0;else print$0}' s1.fin.vcf >s1.title.fin.vcf

6. consensus 序列提取

 ~/bin.bcftoos2/usr/local/bin/vcfutils.pl vcf2fq s3.var.vcf >s3.consensus
bgzip -c s1_var.vcf >s1_var.vcf.gz
~/bin.bcftoos2/usr/local/bin/bcftools index s1_var.vcf.gz
cat ref.fa|~/bin.bcftoos2/usr/local/bin/bcftools consensus s1_var.vcf.gz >consensus/s1.con