Demo entry 6636601

1

   

Submitted by wq on Aug 26, 2017 at 23:04
Language: Perl. Code size: 1.1 kB.

# 对比两个比对软件的输出。
# 使用方法: bash compare.sh
# 在出错的地方停止运行。
set -ue

# 数据文件名。
READ1=r1.fq
READ2=r2.fq

# 参考序列文件。两个比对软件需要分别对它进行索引。
REFS=~/refs/852/NC.fa

# 从基因组中生成模拟reads。
# 编辑错误率并重新运行这个pipeline。
wgsim -N 10000 -e 0.1 $REFS $READ1 $READ2 > mutations.txt

# 将mutations文件转成GFF文件。
cat mutations.txt | awk -F '\t' ' BEGIN { OFS="\t"; print "##gff-version 2" } { print $1, "wgsim", "mutation", $2, $2, ".", "+",  ".", "name " $3 "/" $4 }' > mutations.gff

# 我们需要添加一个read组到mapping中。

GROUP='@RG\tID:123\tSM:liver\tLB:library1'

# 运行bwa并创建bwa.sam文件。
bwa mem -R $GROUP $REFS $READ1 $READ2 > bwa.sam

# 运行bowtie2并创建bow.sam文件。
bowtie2 -x $REFS -1 $READ1 -2 $READ2 > bow.sam

# 调整bowtie2# bowtie2 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x $REFS -1 $READ1 -2 $READ2 > bow.sam
# 对每个sam文件,都转换成bam(sam的二进制)格式。
for name in *.sam;
do
samtools view -Sb $name > tmp.bam
samtools sort -f tmp.bam $name.bam
done

# 删除中间文件。
rm -f tmp.sam tmp.bam

# 修改名字,使得他们不再含有两个后缀。
mv bwa.sam.bam bwa.bam
mv bow.sam.bam bow.bam

# 对数据进行索引。
for name in *.bam
do
samtools index $name
done

# 删除这个程序生成的所有文件。
# rm -f *.bam *.bai *fq *.txt *.sam *.gff

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).