sam和bam格式文件的shell小练习
# 创建名为rna的小环境
conda create -n rna python=2 -c anaconda
# 激活小环境
conda activate rna
# 安装软件
conda install -y sra-tools
# 批量安装
conda install -y sra-tools fastqc trim-galore hisat2 subread multiqc samtools salmon
# 退出小环境
conda deactivate
mkdir -p ~/biosoft
cd ~/biosoft
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
unzip bowtie2-2.3.4.3-linux-x86_64.zip
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq > tmp.sam
# samtools view -bS tmp.sam >tmp.bam
- 统计共多少条reads(
pair-end reads这里算一条)参与了比对参考基因组
把头文件以@开头的去掉，grep -v就是不要，最后20000条，但是左右两条序列算一条，是10000条
- 统计共有多少种比对的类型(即第二列数值有多少种)及其分布。
3）筛选出比对失败的reads，看看序列特征。
cat tmp.sam |grep -v '^@'|less -SN
cat tmp.sam |grep -v '^@'|awk '{if($6=="*")print $10 }'|less -S
- 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID
cat tmp.sam |grep -v '^@'|awk '{if($6=="*")print $1 }'|sort|uniq -c
cat tmp.sam |grep -v '^@'|awk '{if($6=="*")print $1 }'|sort|uniq -c|grep -1
grep -wThe test is that the matching substring must either be at the beginning of the line, or preceded by a non-word constituent character. Similarly, it must be either at the end of the line or followed by a non-word constituent character. Word-constituent characters are letters, digits, and the underscore.
- 筛选出比对质量值大于30的情况（看第5列）
cat tmp.sam |grep -v '^@'|awk '{if($5>30)print $0}'|less -S
- 筛选出比对成功，但是并不是完全匹配的序列
比对成功的就是不带"*"号键的。
cat tmp.sam |grep -v '^@'|awk '{if($6!="*")print $6 }'|less -S
除了匹配的M还包含其他字符串的，比如I、D。
cat tmp.sam |grep -v '^@'|awk '{if($6!="*")print $6 }'|grep "[IDNXSHP]"
- 筛选出inset size长度大于1250bp的
pair-end reads
- 统计参考基因组上面各条染色体的成功比对reads数量
- 筛选出原始fq序列里面有N的比对情况
N在第10列，匹配到第10列含有N的序列,匹配~
awk '{if($10 ~"N")print}'tmp.sam |less -S
- 问题是为什么不用cut呢？因为想看的是整个sam ，如果要是单纯想看含有N的序列的话在fq里就行了呀，但是现在想看的是整个sam 情况
- 另外可以看倒数第2行，虽然有N，但是还是有98个match，然后得到下一题的出题
- 筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况
less -SN tmp.sam|grep -v '^@' |awk '{if($10~N)print}'|awk '{if($6!~"[IDNSHP]")print}'|awk '{if($6!~"*") print}'|less -SN
- sam文件里面的头文件行数
grep '^@' tmp.sam |wc
- sam文件里每一行的tags个数一样吗
cut -f 12-1000 tmp.sam |less -S
不一样
- sam文件里每一行的tags个数分别是多少个
- sam文件里记录的参考基因组染色体长度分别是？
在头文件中可以看到，这个测试sam里只有一条染色体，48000多个字符
- 找到比对情况有insertion情况的
awk '{if($6~"I")print}' tmp.sam |less -S
- 找到比对情况有deletion情况的
awk '{if($6~"D")print}' tmp.sam |less -S
17）取出位于参考基因组某区域的比对记录，比如 5013到50130 区域
less tmp.sam | grep -v '^@'|awk '{if($4>5013 && $4 <50130)print}'|less -S
- 把sam文件按照染色体以及起始坐标排序
awk '{if($4>5013 && $4 <50130)print}' tmp.sam|sort -k4,4|less -S
- 找到 102M3D11M 的比对情况，计算其reads片段长度。
grep 102M3D11M tmp.sam|cut -f 10|wc
- 安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍。
- samtools
