вытяните строку из файла vcf с помощью awk

вытяните строку из файла vcf с помощью awk

04.12.2019 03:05:49 Просмотров 28 Источник

я использую следующий код для обработки числовых данных в таблице vcf.

 cat inputfile | while read row; do
                echo $row > tmp
                originalProb= `awk '{print $1}' tmp`
                probabilityHom1=`awk '{print $2}' tmp`
                probabilityHom2=`awk '{print $4}' tmp`
                numCols=`awk '{print NF}' tmp`

                if [ $numCols -gt 4 ]; then
                        echo "${originalProb}" >> currentRowGenotypes
                elif [ "$probabilityHom1" -gt "$probabilityHom2" ]; then
                        echo "1/1" >> currentRowGenotypes
                elif [ "$probabilityHom1" -lt "$probabilityHom2" ]; then
                        echo "0/0" >> currentRowGenotypes
                elif [ "$probabilityHom1" -eq "$probabilityHom2" ] && [ "$probabilityHom1" -eq 0 ]; then
                        echo "${originalProb}" >> currentRowGenotypes
                else                    
                        echo "het" >> currentRowGenotypes
                fi

        done

        cat tmpHeaders currentRowGenotypes > currentFullCol

входной файл выглядит следующим образом 

1/1     255     231     0
0/1     255     0       152
0/1     255     0       82
0/1     255     0       151
0/1     239     0       31
0/1     255     0       255

по какой-то причине команда awk не распознает первую колонку. есть предложения ?

bash unix bioinformatics
У вопроса есть решение - Посмотреть?

https://stackoverflow.com/questions/59167174/pull-string-from-a-vcf-file-using-awk#comment104557918_59167174
https://stackoverflow.com/questions/59167174/pull-string-from-a-vcf-file-using-awk#comment104558890_59167174
как и выше, но все это может быть одним скриптом awk. Видишь grymoire.com/Unix/Awk.html . Удачи.

Ответы - вытяните строку из файла vcf с помощью awk / pull string from a vcf file using awk

Является ответом!
tshiono

04.12.2019 11:02:14

Это не очень хорошая идея, чтобы создать временный файл только для того, чтобы сделать awk разделите строку на столбцы, потому что:

  • Это приводит к накладным расходам на создание временного файла строка за строкой.
  • Он порождает дочерние процессы несколько раз, чтобы вызвать awk.
  • Причиной ошибки может быть различие синтаксиса между bashи awk.

Вы можете сделать это без использования awk. Пожалуйста, попробуйте следующее:

while read -ra row; do
    originalProb="${row[0]}"
    probabilityHom1="${row[1]}"
    probabilityHom2="${row[3]}"
    numCols="${#row}"

    if (( numCols > 4 )); then
        echo "$originalProb" >> currentRowGenotypes
    elif (( probabilityHom1 > probabilityHom2 )); then
        echo "1/1" >> currentRowGenotypes
    elif (( probabilityHom1 < probabilityHom2 )); then
        echo "0/0" >> currentRowGenotypes
    elif (( probabilityHom1 == probabilityHom2 &&  probabilityHom1 == 0 )); then
        echo "$originalProb" >> currentRowGenotypes
    else
        echo "het" >> currentRowGenotypes
    fi
done < inputfile

cat tmpHeaders currentRowGenotypes > currentFullCol

Как неоднократно предлагали другие, лучшим способом будет писать с awk:

awk '{
    originalProb = $1
    probabilityHom1 = $2
    probabilityHom2 = $4
    numCols = NF

    if ( numCols > 4 )
        print originalProb >> "currentRowGenotypes"
    else if ( probabilityHom1 > probabilityHom2 )
        print "1/1" >> "currentRowGenotypes"
    else if ( probabilityHom1 < probabilityHom2 )
        print "0/0" >> "currentRowGenotypes"
    else if ( probabilityHom1 == probabilityHom2 && probabilityHom1 == 0 )
        print originalProb >> "currentRowGenotypes"
    else
        print "het" >> "currentRowGenotypes"
}' inputfile

cat tmpHeaders currentRowGenotypes > currentFullCol

Надеюсь, это поможет.

Источник
metageni

20.01.2020 04:50:16

Почему бы не использовать Pysam ? Он отлично подходит для разбора BCF/VCF.

Источник

Другие вопросы

Ответы на ваши вопросы

Закрыть X